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

Philippe Collin, David Malec, Yvan Lefevre

LAPLACE, Université de Toulouse, CNRS, INPT, UPS, France.

DOI: 10.4236/aast.2021.61003                

Abstract

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.

Keywords

Partial Discharges, Dielectric Constant, Electric Motor, Finite Elements, Paschen’s Criterion, Ballistic Method

1. 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₂ 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₂ 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

Figure 1. Increasing on board electrical equipment demand in commercial aviation.

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.

2. 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:

$curlE=-\left(\partial B\right)/\partial t$ (1)

$divD=\rho$ (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:

$curlE=0$ (3)

$divD=0$ (4)

Because $E$ is curl free, it can be expressed as a gradient of a scalar potential. That is the electric potential V:

$E=-gradV$ (5)

The electric flux density $D$ is derived from the electric field intensity $E$ and the medium properties:

$D={\epsilon }_0\ast {\epsilon }_r\ast E$ (6)

With ε₀ 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:

${\nabla }^{2}V=0$ (7)

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.

3. 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.:

3.1. 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 on a field line can be calculated step by step with the following equations:

$l\times E=0$ (8)

$l_x\ast E_y-l_y\ast E_x=0$

With l_x and l_y respectively the projections of $l$ along $x$ and $y$ axis in a 2D-system.

From (8) it comes:

$\frac{E_y}{E_x}=\frac{l_y}{l_x}$ (9)

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:

${\Delta }_y={\Delta }_x\ast \frac{E_y}{E_x}$ (10)

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):

$\left(y_{i+1}-y_i\right)=\left(x_{i+1}-x_i\right)\ast \frac{E_{y,i}}{E_{x,i}}$ (11)

Let us define a constant step increase ∆_step at each iteration, so that:

${\Delta }_{step}=\frac{\left(y_{i+1}-y_i\right)}{E_{y,i}}=\frac{\left(x_{i+1}-x_i\right)}{E_{x,i}}$ (12)

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:

$\begin{array}{l}{x}_{i+1}=x_i+{\Delta }_{step}\ast E_{x,i}\\ y_{i+1}=y_i+{\Delta }_{step}\ast E_{y,i}\end{array}$ (13)

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₀, y₀) of a field line. Then, the next point on the line M₁(x₁, y₁) is computed by iteration from point Q. The iteration is backward because the starting points are on the end of the field lines:

$\begin{array}{l}{x}_1=x_0-{\Delta }_{step}\ast E_{x,0}\\ y_1=y_0-{\Delta }_{step}\ast E_{y,0}\end{array}$ (14)

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 3 displays the field lines after one iteration backward from starting points. It can be seen that the M₁ points do not coincide with a mesh node. The electric field orientation vectors are interpolated on the extra M₁ points using the interp 2 Matlab function [11]. 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:

$\begin{array}{l}{x}_{i+1}=x_i-{\Delta }_{step}\ast E_{x,i}\\ y_{i+1}=y_i-{\Delta }_{step}\ast E_{y,i}\end{array}$ (15)

The computed field lines after 11 iterations are displays on Figure 3.

Figure 2. Left: System of two infinite cylinders oppositely charges, field of study represents one quarter of the system, mesh (dots) and half conductor contour (blue), right: Electric field orientation vectors and starting points (red circles).

Figure 3. Left: Field lines after a first iteration backward, right: Computed field lines after 11 iterations backward.

3.2. 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:

$\delta \varphi ={\displaystyle {\int }_a^{b}E\cdot n\cdot \text{d}s}$ (16)

where n is the normal vector and s is the path between the lines a and b. This path is chosen as the contour of scalar potential V. Such as in Whittaker’s method, a uniform mesh is applied. The electric field is calculated on the whole meshed domain using linear interpolations. In a simple case, all the field lines pass through a single V contour. The contour results in a series of points $s_i=\left\{x_i,y_i\right\};i=1,\cdots ,I$. The total flux through the contour of potential V is:

${\varphi }_V={\displaystyle {\int }_{s_1}^{s_I}E\text{d}s}=\frac{1}{2}{\displaystyle {\sum }_{i=1}^I\left(E_i+E_{i+1}\right)\left|s_{i+1}-s_i\right|}$ (17)

From (17) to (21) E refers to the normal electric field component. Figure 4 displays the contour V discretization in a uniform meshed domain. Here the total number of points is I = 5. The field lines starting points on contour V results in a second series of points $s_j=\left\{x_j,y_j\right\};j=1,\cdots ,N$. Such points satisfy the following equation:

${\displaystyle {\int }_{s_j}^{s_{j+1}}E\text{d}s}=\delta \Phi$ (18)

Starting from the s₁ point of the contour V, the s_i contour point which is in proximity of the next field line starting point is identified by this equation:

${\displaystyle {\int }_{s_1}^{s_i}E\text{d}s}<\delta \varphi <{\displaystyle {\int }_{s_1}^{s_{i+1}}E\text{d}s}$ (19)

Figure 4. Contour discretization (red segments) in a uniform mesh (black).

The s_j field line starting point on contour V segment s_i, s_i+1 is defined as:

${\displaystyle {\int }_{s_i}^{s_j}E\text{d}s}=\delta \varphi -{\displaystyle {\int }_{s_1}^{s_i}E\text{d}s}={\varphi }_0$ (20)

Equation (20) can be expressed as:

$\frac{1}{2}\left(E_i+E_j\right)\left|s_j-s_i\right|={\varphi }_0$ (21)

Let $a=\frac{\left|s_j-s_i\right|}{\left|s_{i+1}-s_i\right|}$ be the fractional distance from contour point s_i to field line starting point s_j. As the electric field $E$ is linear along the grid, one finally gets (22) [16]:

$s_j=\left\{x_j,y_j\right\}$ (22)

$s_j=\left\{\left(1-a\right)\ast x_i+a\ast x_{i+1},\left(1-a\right)\ast y_i+a\ast y_{i+1}\right\}$ (23)

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₁ and V₂ respectively. A field line starting from V₁ contour is identified by $l_{n1}^1,n1=1,\cdots ,N1$. The same procedure described in the simple case is applied to V₁ contour. However, the intersections of $l_{n1}^1$ field lines with V₂ contour have to be tracked. To do that, each segment of $l_{n1}^1$ is checked to verify whether or not it intersects with V₂ contour. The intersection points are added to V₂ contour points series. Now, one can apply the same procedure with V₂ contour adding some steps:

• Step 1: choose as a starting point of field lines from V₂ contour ( $l_{n2}^2,$ $n2=1,\cdots ,N2$) one of the intersection point of a $l_{n1}^1$ field line with V₂ contour;

• Step 2: integrate from that point until (a) another intersection point of a $l_{n1}^1$ field line with V₂ contour is reached or (b) the integral exceeds the fixed flux quantity δϕ. In case (b), the next point as to be determined the same way s_j point is in the simple case;

• Step 3: repeat step 2 until the whole V₂ 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₁ and V₂ using the presented algorithm.

4. 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.

4.1. 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.

$\{\begin{array}{l}{x}_e\left(u,v\right)={\displaystyle {\sum }_{i=1}^8{N}_i\left(u,v\right)\cdot x_e\left(i\right)}\\ y_e\left(u,v\right)={\displaystyle {\sum }_{i=1}^8{N}_i\left(u,v\right)\cdot y_e\left(i\right)}\\ V_e\left(u,v\right)={\displaystyle {\sum }_{i=1}^8{N}_i\left(u,v\right)\cdot V_e\left(i\right)}\end{array}$ (24)

With:

Figure 5. Electric field lines from two infinitely long wires of opposite charge in free space. V₁ and V₂ are equipotential lines. The solid field lines were drawn from V₁. The dashed field lines were drawn from V₂.

$\begin{array}{l}{N}_i\left(u,v\right)=1\text{at}\text{node}i\\ N_i\left(u,v\right)=0\text{elsewhere}\end{array}$

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):

$\left(\begin{array}{c}\frac{\partial }{\partial u}\\ \frac{\partial }{\partial v}\end{array}\right)=\left[\begin{array}{cc}\frac{\partial x}{\partial u}& \frac{\partial y}{\partial u}\\ \frac{\partial x}{\partial v}& \frac{\partial y}{\partial v}\end{array}\right]\left(\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\end{array}\right)=J\left(\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\end{array}\right)$ (25)

With:

$J=\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]$

With J the Jacobian matrix. It is computed from the known shape functions partial derivatives in the local system when putting (24) in (25). Then, the partial derivatives of a function in the global system (x, y) can be expressed from its partial derivatives in the local system (u, v):

$\left(\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\end{array}\right)=\left[\begin{array}{cc}\frac{\partial u}{\partial x}& \frac{\partial v}{\partial x}\\ \frac{\partial u}{\partial y}& \frac{\partial v}{\partial y}\end{array}\right]\left(\begin{array}{c}\frac{\partial }{\partial u}\\ \frac{\partial }{\partial v}\end{array}\right)=I\left(\begin{array}{c}\frac{\partial }{\partial u}\\ \frac{\partial }{\partial v}\end{array}\right)$ (26)

With:

$I=\left[\begin{array}{cc}{I}_{11}& I_{12}\\ I_{21}& I_{22}\end{array}\right]$

$I=J^{-1}$

It is matrix I which is used in practice because the data have to be expressed in the global system. It is computed as follow:

$I=J^{-1}=\frac{1}{\det \left(J\right)}\left[\begin{array}{cc}{J}_{22}& -J_{12}\\ -J_{21}& J_{11}\end{array}\right]$ (27)

With:

$\det \left(J\right)=J_{11}{J}_{22}-J_{12}{J}_{21}$

The electric field components on the nodes of an element e can finally be expressed in the global system by combining (24) and (26):

$\{\begin{array}{l}{E}_{x,e}=-\frac{\partial V_e}{\partial x}=-{\displaystyle {\sum }_{i=1}^8\left(I_{11}\frac{\partial N_i}{\partial u}+I_{12}\frac{\partial N_i}{\partial v}\right)\cdot V_e\left(i\right)}\\ E_{y,e}=-\frac{\partial V_e}{\partial y}=-{\displaystyle {\sum }_{i=1}^8\left(I_{21}\frac{\partial N_i}{\partial u}+I_{22}\frac{\partial N_i}{\partial v}\right)\cdot V_e\left(i\right)}\end{array}$ (28)

The electric flux components are computed from electric field components using (6).

4.2. 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:

$z\left(t\right)=a_{1,z}+a_{2,z}\ast t+a_{3,z}\ast t^2$ (29)

With:

$t\in \left[0,1\right]$

Figure 6 displays (29) along a parameterized edge. For instance, on Figure 6 the interpolation is done on the edge containing the nodes (1, 2, 3). The nodes are parameterized according to the local numbering order: t = 0 for node 1, t = 0.5 for node 2 and t = 1 for node 3. By solving (29) in term of voltage it is possible to derive the polynomial coefficients:

Figure 6. Left: Equipotential line V intersecting an element edges, right: Illustration of a second order polynomial interpolation over a parameterized edge.

$V\left(t\right)=a_{1,V}+a_{2,V}\ast t+a_{3,V}\ast t^2$ (30)

With:

$\{\begin{array}{l}V\left(0\right)=a_{1,V}=V_0\\ V\left(0.5\right)=a_{1,V}+a_{2,V}\ast 0.5+a_{3,V}\ast 0.25=V_{0.5}\\ V\left(1\right)=a_{1,V}+a_{2,V}+a_{3,V}=V_1\end{array}$

In the considered example, (V₀, V_0.5, V₁) 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.

$V\left(t_0\right)=V=a_{1,V}+a_{2,V}\ast t_0+a_{3,V}\ast t_0^2$ (31)

With V the voltage of the considered equipotential line V. The obtained parameter t₀ 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.

4.3. 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 illustrates the barycenter reference attached to an equipotential line V. Point G is the barycenter. A point M located on the equipotential line is identified by its curvilinear abscissa s(M). A flux function ϕ(M) of points on each equipotential is defined as the flux per meter crossing the line between the starting point Q and point M:

$\varphi \left(M\right)={\displaystyle {\int }_Q^M{D}_n\text{d}s}$ (32)

Figure 7. Left: Barycenter system associated to an equipotential line V, right: Flux tubes computation.

D_n is the normal electric flux density on the equipotential line, ds is the elemental curvilinear length of the equipotential line. Each point M on the line is parameterized by the curvilinear abscissa s(M) and:

$\{\begin{array}{l}s\left(M\right)={\displaystyle {\int }_Q^M\text{d}s}\\ x\left(M\right)=x\left(s\right)\\ y\left(M\right)=y\left(s\right)\\ \varphi \left(M\right)=\varphi \left(s\right)\end{array}$ (33)

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 δϕ:

$\delta \varphi ={\displaystyle {\int }_{P_i}^{P_{i+1}}{D}_n\text{d}}s$ (34)

The properties of these points are:

$\{\begin{array}{l}{P}_1=Q\\ \varphi \left(P_1\right)=0\\ \varphi \left(P_i\right)=\varphi \left(P_i-1\right)+\delta \varphi =\left(i-1\right)\delta \varphi \end{array}$ (35)

All the points P_i on the equipotential line are determined by inversing the previous functions:

$s\left(P_i\right)={\varphi }^{-1}\left(\left(i-1\right)\delta \varphi \right)$ (36)

Figure 7 illustrates the flux tubes computation using (33), (34), (35) and (36). The determination of points P_i is done on all equipotential lines. If the starting point Q on each equipotential line is correctly chosen, then the electric flux lines are defined by iso-ϕ lines. The starting point has to be chosen on a field line crossing all the equipotential lines present in the domain.

5. 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 stator slot in the presence of a default (air bubbles or bad impregnation) in the surrounding impregnation resin. The copper areas are obviously not modelled since there is no electric field in conductor materials. Boundary conditions (i.e.: voltages) are applied on the copper areas contour. The air and the enamel are considered charge carrier free. The wires are considered infinitely long in the machine active length dimension and the voltage drops are neglected. Due to invariance along this dimension and symmetries, the problem is reduced to a two dimensional (2D) electrostatic problem with only a quarter of the wires being modelled. The described 2D problem is displayed on Figure 8. In such electrostatic problem, the materials are characterized with their dielectric constant. 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.

Figure 8. Enameled magnet wires 2D-electrostatic problem.

Table 1. Wire parameters.

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.

5.1. 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.

5.2. 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

Figure 9. Finite elements mesh, white: copper, purple: enamel, blue: air.

Figure 10. Total electric fluxes crossing each equipotential line.

Figure 11. Equipotential lines (colored dotted lines) and electric field lines (colored solid lines), only the enamel external contour is represented (black solid lines).

implanted on Matlab function stream2 [11]. That function requires a uniform mesh. However, the finite elements mesh on Figure 9 is not uniform. A second mesh is thus created on Matlab using meshgrid function [11]. Axis vectors X_a, Y_a are used as inputs of meshgrid. They are built from X, Y finite elements nodes coordinates lists such that:

• $X_a\in \left[X_{\min },X_{\max }\right]$ and $Y_a\in \left[Y_{\min },Y_{\max }\right]$

• 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.

5.3. 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⁻⁴. 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

Figure 12. Electric field lines in air computed with Matlab ballistic method (solid black lines), colored solid lines represent the external enamel contours.

Figure 13. Electric field amplitude [V/m], applied voltage 1000 V peak.

Figure 14. Computed voltage drops along field lines in air (p = 1 bar) obtained with the two compared methods, γ = 9 × 10⁻⁴.

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.

6. 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.

Acknowledgements

This project has received funding from the Clean Sky 2 Joint Undertaking under the European Union’s Horizon 2020 research and innovation program under grant agreement No 715483.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	(1992) Agenda 21 (Rio Declaration on Environment and Development). 79.
[2]	International Energy Agency (2018) IEA World Energy Outlook 2018.
[3]	International Transport Forum, OECD iLibrary (2017) ITF Transport Outlook 2017.
[4]	Organization of the Petroleum Exporting Countries (2018) OPEC World Oil Outlook 2040.
[5]	World Health Organization (2019) WHO Health and Sustainable Development: Transport and Health Risks.
[6]	Rodrigue, J.P., Comtois, C. and Slack, B. (2017) The Geography of Transports Systems. 4th Edition, Routledge, New York. https://doi.org/10.4324/9781315618159
[7]	Robinson, J. (2019) Collins Aerospace Building $67M Electric Aircraft Lab. Wings Mag.
[8]	Sarlioglu, B. and Morris, C.T. (2015) More Electric Aircraft: Review, Challenges, and Opportunities for Commercial Transport Aircraft. IEEE Transactions on Transportation Electrification, 1, 54-64. https://doi.org/10.1109/TTE.2015.2426499
[9]	Hemmati, R., Wu, F. and El-Refaie, A. (2019) Survey of Insulation Systems in Electrical Machines. 2019 IEEE International Electric Machines & Drives Conference (IEMDC), San Diego, May 2019, 2069-2076. https://doi.org/10.1109/IEMDC.2019.8785099
[10]	Szczepanski, M. (2019) Development of Methods Allowing the Test and the Comparison of Low-Voltage Motors Insulation Systems Running under Partial Discharges (Fed by Inverter). 214.
[11]	MATLAB, 9.3.0.713579 (R2017b). The MathWorks Inc.
[12]	Ren, Z. (2003) 2-D Dual Finite-Element Formulations for the Fast Extraction of Circuit Parameters in VLSI. IEEE Transactions on Magnetics, 39, 1590-1593. https://doi.org/10.1109/TMAG.2003.810548
[13]	Parent, G., Duchesne, S. and Dular, P. (2017) Determination of Flux Tube Portions by Adjunction of Electric or Magnetic Multivalued Equipotential Lines. IEEE Transactions on Transportation Electrification, 53, 1-4. https://doi.org/10.1109/TMAG.2017.2661401
[14]	Parent, G., Rossi, M., Duchesne, S. and Dular, P. (2019) Determination of Partial Discharge Inception Voltage and Location of Partial Discharges by Means of Paschen’s Theory and FEM. IEEE Transactions on Transportation Electrification, 55, 1-4. https://doi.org/10.1109/TMAG.2019.2902374
[15]	Whittaker, T.M. (1977) Mon Weather.
[16]	Horowitz, E.J. (1990) A General Field-Line Plotting Algorithm. Journal of Computational Physics, 4, 418. https://doi.org/10.1063/1.168384
[17]	Ansys® Academic Research Mechanical, Release 16.1.
[18]	Dhatt, G. and Touzot, G. (1981) Une présentation de la méthode des éléments finis, Les Presses de l’Université Laval. Québec, et Maloine, Paris.
[19]	Collin, P., Malec, D. and Lefevre, Y. (2019) About the Relevance of Using Paschen’s Criterion for Partial Discharges Inception Voltage Estimation When Designing the Electrical Insulation System of Inverter Fed Motors. Electrical Insulation Conference (EIC), Calgary, 16-19 June 2019, 513-516. https://doi.org/10.1109/EIC43217.2019.9046558