Application of rotating coupling time-step ﬁnite element method in synchronous generators’ internal faults simulation

A precise simulation of the internal faults of synchronous generators is very important for designing the main protection scheme; therefore, it is necessary to propose an accurate model. The time-step ﬁnite element method is quite qualiﬁed for analysing problems related to the transient electromagnetic ﬁeld, including the internal faults of synchronous generators. This paper presents a new time-step method, which employs a rotation coupling technique to handle rotor motion. By using this technique, the rotor remains stationary during the process of the simulation, and the motion problem becomes simple and can be conveniently implemented with a programming process. Simulation and experimental results for the internal faults of an experimental machine are compared to verify the accuracy of the rotation coupling time-step method. Moreover, a 300 MW salient pole synchronous generator is used as an example, and the transient process of this generator during internal faults is simulated. The air gap ﬂux density, fault currents, damper windings eddy current losses, and core dynamic electromagnetic force are investigated in detail to reveal the fault characteristics more clearly and provide a basis for designing


INTRODUCTION
The internal faults of synchronous generators can result in serious damage to generators and even the power system [1]. Hence, it is very important to design an adequate protection scheme for synchronous generators. However, the design of an adequate protection scheme is dependent on the accurate calculation of the transient fault currents [2]. Many methods have been employed to analyse the internal faults of synchronous machines, such as two-reaction theory [3], the direct phase representation [4], the winding function theory [5,6], and the multi-loop circuit method [7,8]. With the improvements in computer performance, the application of the time-step finite element method has increasingly widened. This method can accurately simulate the transient process of internal faults because it can account for all of the space harmonics, material non-linearities, and the influence of slotting [9][10][11][12][13]. However, problems with the rotor motion are inevitable when using the time-step finite element method.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. A number of methods have been proposed to deal with rotor motion problems. The moving boundary technique is simple and efficient, which is widely used in the transient simulation analysis of machines [14,15]. The moving band technique is also widely used in motion problems [16]. However, the re-meshing of moving parts may cause distortions in the air gap elements. The air gap macro-element method introduced by Razek treats the air gap as a single entity to analyse rotating electrical machines [17]. The no mesh method is usually employed when the motion situation is complex [18]. In three-dimensional time-step finite element problems, the lockstep method [19] and Lagrangian-multiplier interpolation method [20] are preferentially used. The former is suitable for a fixed step motion, and the latter is competent for handling arbitrary step problems.
There are thousands of possible short-circuit faults in the stator windings of large synchronous generators [21]. However, the design of the main protection scheme of generators needs to be accurately simulated for each short-circuit fault, which is very important to choose adequate relay type and trip level [22]. Obviously, the amount of calculation and work are very huge. Therefore, it is very important to save simulation time and implement a programming process. This paper proposes the rotation coupling technique first. By coupling the magnetic vector potential of nodes located in the sliding interface, the continuity of the stator and rotor magnetic field can be maintained at each step. The rotor can remain stationary during the simulation. This technique simplifies the motion problem, and the time-step finite element problem can be conveniently implemented with a programming process. Stator winding currents obtained from simulated results are compared with experimental results to verify the accuracy of the rotation coupling technique. Internal fault problems for a 300 MW salient pole synchronous generator are presented as an example, and the air gap flux density, fault currents, damper bars eddy current losses, and core dynamic electromagnetic force are simulated using the rotation coupling time-step finite element method, thereby revealing the fault characteristics more clearly and providing a basis for designing protection schemes.

Magnetic-field equations
The following assumptions are made: 1. For low-frequency electromagnetic field problems, the ratio of the displacement current to the conduction current is less than 10 −7 , so the influence of the displacement current is ignored [23]; 2. because the internal short-circuit simulation focuses on the transient process in the first few cycles after the fault occurs, at this time, the temperature of the winding caused by the increase of the fault current does not change much, so the influence of the temperature on the conductivity of the winding is ignored; and 3. the core is considered to be infinitely long in the z-direction [23], and the end winding resistance and reactance are included in the external circuits; therefore, the magnetic field of the machines is considered to be a two-dimensional distribution.
It is generally considered that the electromagnetic field in the electric machines is a quasi-steady field, and the electromagnetic relationship can be obtained by the Maxwell equation where H is the magnetic field strength, J is the current density in the stator coils, excitation coils and damping bars, A is the vector magnetic potential, and σ is the conductivity. Ignoring the hysteresis effect of ferromagnetic materials, the relationship between H and A is where ν is the reluctivity. Assumption (3) shows that the vector magnetic potential and current density have only axial components, that is, A = A(x,y)⋅z, J = J⋅z. Substituting Equations (2) into (1), the transient electromagnetic field equation can be obtained as where Γ 1 and Γ 2 are the outer and inner boundaries, respectively. Control Equation (3) is discretised using the Galerkin method and a matrix equation is obtained as where K and M are the coefficient matrices associated with the node coordinates, p is the differential operator, I b is the coil currents matrices of the stator, field, and damper coils, respectively, and C is the correlation matrix between the coil currents and nodes.

New method for handling the rotor motion
The rotor motion is simulated by the rotation coupling technique proposed in this paper. In order to implement this technique, the air gap is divided into at least two layers. The part of the grid in the air gap region is shown in Figure 1. There are two rotation lines in the middle of air gap, which belong to the upper and lower gap elements. The nodes of the rotation lines To simulate the motion of the rotor, the magnetic vector potential of the node located in the rotation lines is coupled in a certain order and in turn. For the starting moment of the rotor rotation, that is, t 0 = 0, the nodes S(1) and R(1) on the two rotation lines are coincident, nodes S(2) and R(2) are coincident and so on. In this case, as shown in Figure 2(a), the magnetic vector potential of the coincident nodes should be coupled as follows: where N is the total number of nodes on one rotation line.
When the simulation time is t 1 = t 0 +Δt, where Δt is the simulation step, assuming that the rotor actually rotates θ angle, the angular speed ω of the generator is = ∕Δt (6) The constraint condition between N, ω and is where P is the number of pole pairs. To simulate the motion of the rotor, at time t 1 , the magnetic vector potential of the nodes should be decoupled first, which Using the rotation coupling technique, the rotor model is maintained stationary during the simulation, so the time of handling rotor motion is saved and the pre-process is simplified. In addition, because the node labels are stored in arrays, the rotation coupling technique can be implemented with a programming process.

Time-step finite element equations
The internal fault problem of synchronous generators is taken as an example. Figure 3 shows a schematic of the external circuit-coupled finite element model. R s and L s represent the resistance and leakage inductance of the stator end windings, respectively. The effective length of the Finite element (FE) region is equal to the core length, and the slot resistances and leakage inductances of the windings are included in this region. R t and L t represent the equivalent resistance and inductance of the transformer, respectively. According to the external circuit model in Figure 3, the following circuit equation can be obtained where U, Ψ, and I are the matrices of the loop voltages, loop flux linkages, and loop currents, respectively, and R is the resistance matrix.
The number of loops in case of the internal fault of generators is (10) where N s is the number of stator loops and N s = 3a; a is the number of parallel branches; N f is the number of field loop and N f = 1; N d is the number of damper loops and N d = 2P⋅N c ; N c is the number of damping bars per pole.
The relationship between the loop current matrix and the coil current matrix is where G is the incidence matrix of the loop current and branch current, and the calculation method of G is with reference to [22]. The flux linkage Ψ is composed of the flux linkage in the axial part of the generator and the end leakage flux linkage where Ψ L is the flux linkage in the axial part of the generator, Ψ E is the end leakage flux linkage. The flux linkage in the axial part of the generator can be calculated by the finite element equation, and the Ψ L can be written as where l ef is the effective length of the generator. The end leakage flux linkage Ψ E can be written as where M E is end leakage inductance, and the calculation method of M E is with reference to [22]. Substituting Equations (12) into (9) results in (15), the control equation can be written as follows: The motion equation of the rotor can be expressed as follows: where T e is the electromagnetic torque, T L is the load torque, J is the moment of inertia, and Ω is the rotor speed. The Equations (16) and (17) constitute the fundamental equations of the step-time finite element method.

EXPERIMENTAL VERIFICATION
In order to verify the accuracy of the rotating coupling timestep finite element method, the three-phase sudden short-circuit fault and the stator winding internal fault in a 30 kVA synchronous generator are analysed. The time-step finite element model consists of 103,904 nodes and 51,027 triangular elements. The laboratory testbed is shown in Figure 4. The synchronous generator is connected to an infinite bus through a transformer.
The rated voltage and current of this generator are 400 V and 43.3 A, respectively. This generator has six poles and two parallel branches for each phase. The measured and simulated waveforms of the stator phase currents I A , I B , and I C under the three-phase short-circuit fault condition are compared as shown in Figures 5-7. In order to limit the short-circuit current, the excitation current I f = 0.178 A (rated excitation current I fN = 0.958 A), and the generator operates in a no-load state.
The internal fault occurs in the end region between coils no. 2 and no. 3 of branch A1, and the generator operates at half load. The fault currents in branches A1 (I A1 ), B1 (I B1 ), and C1 (I C1 ) were simulated and measured by current transformers. Figures 8-10 show the comparisons between the measured and simulated current waveforms of I A1 , I B1 , and I C1 . The internal fault occurred in 0.02 s. It can be seen that the branch currents are basically symmetrical before the fault, and the amplitude of the branch current is 15.3 A at half load. However, the branch currents are no longer symmetrical after the fault. Because only one coil is short-circuited, the amplitude of each branch current does not significantly change.
In summary, for either normal or fault (three-phase sudden short-circuit fault and internal fault) operation, the simulated waveforms agree well with measured waveforms; thus, the accuracy of the rotating coupling time-step finite element method is verified.

EXAMPLE OF INTERNAL-FAULT SIMULATION
The rotating coupling time-step finite element method is used to simulate the internal faults of a 300 MW salient pole synchronous generator. The major parameters of this generator are listed in Table 1. Figure 11 shows the models of the synchronous generator for internal fault simulations. The FE model consists of 121,679 nodes and 61,182 triangular elements. To ensure calculation accuracy, the element sizes of the air gap and conductors are suitably controlled. For example, the sizes of the air gap elements are between 8-10 mm. The magnetic permeability of the stator and rotor core is non-linear. The relative permeability of the stator, field, and damper wind- ings is equal to 1, and the electrical resistivity is equal to 1.75 × 10 -8 Ω⋅m. The Newton-Raphson method is used to solve the non-linear equations [24]. The simulation time-step size is 0.4 ms.
The internal fault occurs in the end region between coils no. 11 and no. 12 of branch A1, so the fault is the inter-turn shortcircuit and only one coil is short-circuited. The fault circuit diagram is with reference to Figure 3. Figure 12 shows the fault transient current waveforms of I A1 and I k . I A1 is the branch current of branch A1, and I k is the short-circuited loop current. As shown in Figure 12, because only one coil is short-circuited, accounting for 0.47% of the total coil number of A1 branch, the amplitude of I A1 is 2.69 kA, which is only 0.25 times the rated current. However, the amplitude of the short-circuited loop current I k is 98.76 kA, which is 9.22 times larger than the rated current, and can therefore result in serious damage to the generator and its associated system.
Assuming that the short-circuit loop current I k contains only the fundamental component, the expression of I k is where I m is the amplitude of I k , and φ is the initial phase angle of I k .
A pulsating magnetomotive force (MMF) is built by the short-circuit loop current as shown in Figure 13. The fundamental harmonic expression for the pulsating MMF f 1 is  where k is the number of coils that are included in the shortcircuit loop, N k is the coil turns, k w1 is the winding factor, θ is the electrical angle with a range, 0-360 0 , from point A to point B.
The pulsating MMF distorts the air gap magnetic field, resulting in the damper winding eddy current loss and core dynamic electromagnetic force sharp increase. Therefore, we must make an in-depth study of the fault characteristics caused by the internal fault, providing a basis for the optimal design of the generators core and windings structure and the configuration of the main protection scheme. Figure 14 shows the air gap flux density of the entire generator range under the internal fault condition. Because only one coil is short-circuited, the pulsating magnetic field covers 360 • electrical angles, so the air gap magnetic field of only a pair of poles is severely distorted. Figure 14(b) presents the harmonic content of the air gap flux density. The 35th and 37th harmonics are first-order stator tooth harmonics. Obviously, the air gap flux density of the no. 1 pair of poles is severely distorted and the amplitudes of the fundamental harmonics are reduced; that is, at this instant, no. 1 pair of poles are located in the fault area of the stator windings. The amplitude of the fundamental harmonic magnetic field is 0.56 T, and the amplitude of the 3 rd harmonic magnetic field is 0.41 T, accounting for 73% of the fundamental harmonics. Moreover, some even harmonics appear, under the    Figure 15 shows the transient current waveforms of damper bars, and the damper bars labels are with reference to Figure 11. The damper bar currents vary periodically with a period of 0.14 s, which is equal to the time taken for one rotation of the rotor.

Eddy current loss of damper winding
In each period, due to the influence of the short-circuit loop current, when the damper bars pass through the fault area, their currents reach high peak values. The maximum value of the damper bar eddy current is 36.2 kA, which is 3.38 times of the rated stator current. In addition, the currents of the damper bars on both sides of the magnetic pole are significantly greater than where σ is the conductivity of damper bars. The damper bar element eddy current and eddy current loss can be calculated as where L b is the axial length of damper bars, Δ e is the area of element e. Therefore, the eddy current loss of one damper bar is (23) where N is the total number of elements of one damper bar. Figure 16 shows the distribution of element eddy current losses of the damper winding. After the internal fault occurs, the distribution of eddy current losses of damper bars is extremely uneven. The eddy current losses of the no. 1 and no. 9 damper bars on both sides of the no. 5 pole increase sharply; thus, we can conclude that no. 5 pole is located in the fault area, at that time, and the eddy current losses of damper bars in the non-fault area are very small. By calculation, the maximum instantaneous eddy current loss of a single damper bar almost reaches 620 kW, 23,220 times larger than the eddy current loss during no-load operation. Therefore, if the generator is operated in the internal fault state for a long time, it may cause damage to the damper bars on both sides of the pole.

Dynamic electromagnetic force of stator core
Due to the difference in permeability, the electromagnetic force will be produced at the interface between the core and air gap or conductors of the generator, and the direction of the electro- magnetic force is the normal direction of the interface. When an inter-turn short-circuit occurs, because the short-circuit current in the fault coil is very large, it will inevitably lead to an increase in the magnetic field around the fault coil and serious distortion. As a result, the local electromagnetic force of the core around the fault coil increases sharply. Figure 17 shows the distribution of the magnetic flux density of the core around the fault coil. Assume that the upper bar of the fault coil is located in slot 1#, whereas, the lower bar is in slot 19#. As shown in Figure 17, taking the 1# stator tooth as an example, the electromagnetic force on both sides of the tooth wall is calculated. The reference direction of dynamic electromagnetic force and the distribution of key points (P 1 -P 8 ) on the stator tooth are shown in Figure 18. According to the Maxwell stress method, the electromagnetic force on the tooth wall of the stator core can be calculated as where f n and f t are the normal component and tangential component of the core electromagnetic force; μ Fe and μ 0 are the magnetic permeability of the stator core and air gap, respectively; B n and H t are the normal component of the magnetic flux density and the tangential component of the magnetic field strength at the interface between the stator core and air gap, respectively.

FIGURE 19
The change of electromagnetic force of key points and I k with time

FIGURE 20
Electromagnetic force of key points of 1# stator teeth Figure 19 comprehensively shows the changes of the shortcircuit loop current I k and the electromagnetic force of key points with time, where the fault occurred at 0.02 s. As can be seen from Figure 19, before the fault occurred, the amplitudes of the electromagnetic force at the key points P 1 and P 5 are small and are almost equal. The amplitudes of the electromagnetic force at key points P 1 and P 5 are 0.193 × 10 6 N/m 2 and 0.199 × 10 6 N/m 2 , respectively. However, due to the influence of the pulsating magnetic field produced by I k , after the fault, the amplitudes of the electromagnetic force at the key points P 1 and P 5 increase sharply, and the distributions are no longer symmetrical. The electromagnetic force at key point P 1 , next to the fault coil, is much larger than the electromagnetic force at key point P 5 , where the amplitude of the electromagnetic force at point P 1 is 1.04 × 10 6 N/m 2 , which is 5.39 times larger than that before the fault; the amplitude of the electromagnetic force at key point P 5 is 0.545 × 10 6 N/m 2 , 2.74 times larger than that before the fault. Figure 20 shows the spatial distribution of electromagnetic force at key points P 1 -P 4 on the tooth wall after the internal fault. Obviously, the electromagnetic force near the tooth top is much larger than the electromagnetic force near the tooth bottom. Among them, the electromagnetic force amplitude at key points P 1 -P 4 are 1.04 × 10 6 N/m 2 , 0.918 × 10 6 N/m 2 , 0.258 × 10 6 N/m 2 and 0.110 × 10 6 N/m 2 , respectively. Therefore, on the whole, we must pay attention that, after the internal fault, the sharp increase of the electromagnetic force on the stator teeth may cause damage to the generator core. In addition, the maximum value of electromagnetic force is concentrated on the teeth top, which is next to the fault coil.

CONCLUSION
The rotation coupling technique is proposed in this paper to overcome the problems associated with rotor motion. At each time step, instead of rotor rotation, it is only necessary to couple the magnetic vector potential of the nodes located on the rotation lines in a certain order and in turn. Using the rotation coupling technique, the pre-processing of the time-step finite element method is simplified, and the simulation efficiency is improved.
The internal faults of a 300 MW salient pole synchronous generator are simulated using the rotation coupling time-step finite element method, and the following conclusions may be drawn. First, the short-circuit loop current I k is extremely large, which can reach about 10 times larger than the rated current. The pulsating MMF produced by the short-circuit loop current distorts the air gap magnetic field, resulting in damper eddy current losses and dynamic electromagnetic force of core sharp increase. Second, the eddy current losses of the damper bars located in the fault area and on both sides of the magnetic pole are considerably greater than the others. Finally, the maximum value of electromagnetic force is concentrated on the teeth top which is next to the fault coil.