A Comprehensive Analysis on Transient Electromagnetic Force Behavior of Stator Windings in TurboGenerator

This paper presents a comprehensive analysis on the transient electromagnetic force behavior of the stator windings in a QFSN600-2YHG type turbo-generator. Different from other studies, this paper investigates not only the distribution regularities of the resultant force and the force density, but also the force harmonic characteristics, and the mechanical responses which will cause sensitive impact on the insulation wearing. The whole work is generally based on a proposed simplified model and the 3D finite element coupling calculation. The simplified model contains two parts. The first part is the theoretical model that employs the approximate solution of the image current to obtain the analytic formula of the electromagnetic force on the end windings conveniently. The second part is the FEA model that employs only the end windings and one-tenth of the stator core to save the calculating memory and, meanwhile, obtain the qualified electromagnetic force as well as the mechanical response. It is shown that the nose-top, the connection point between the line part and the end part, and the middle of the involute are the three most dangerous positions of the endwinding to sustain serious insulation wearing. Moreover, the winding, which endures themaximum mechanical response, is neither always consistent with the one that has the largest resultant electromagnetic force nor directly in accordance with the winding that affords the most intensive electromagnetic force density. The findings in this paper will be beneficial for the insulation monitoring and the manufacturing improvement on the stator windings.


Introduction
Stator windings are key parts in the generator for energy transformation and electricity output.However, intensive electromagnetic force acting on the windings for a long term may result in degrading the insulation.As the capacity of the generator becomes larger and larger, the electromagnetic forces acting on the stator windings, especially on the end regions, will be increased greatly [1].Therefore, the end winding is easier to suffer wearing and damage than the line part because of the larger forces as well as its cantileverbeam structure [1,2].The excessive forces on the end winding will result in severe vibrations, mal-strains, overstresses, and loose of fasteners [3].After a long-period performance, it is highly probable that the insulation of the winding on some particular parts will be destroyed due to the dynamic vibration friction, overstatic stress wearing, and serious strain degrading [4].
Scholars have carried out many studies on the electromagnetic forces of the end windings.Also, they have gained plenty significant achievements, which founded a good basis for the design and the condition monitoring on the generator windings.The electromagnetic force problem, as well as the significance of the force control, was proposed early by C. A. Bucci et al. in 1971 [5].Hereafter, the finite element analysis (FEA) method and the corresponding case study of the electromagnetic force on the stator end windings under steady state [6] and transient conditions [7] in large turbine generators were presented, respectively.In addition, the integral equation method for calculating the electromagnetic forces under symmetric, asymmetric, sudden short-circuit, and various power factor conditions [8], was proposed as well.It is found that the quasi-three-dimensional (quasi-3D) finite-element method is able to take advantages over the conventional 3D FEA in memory size and time efficiency [9].Significantly, besides the rated load condition, numerical 2 Mathematical Problems in Engineering simulations on the electromagnetic force densities in the three-phase short circuit condition [10] and the leading power factor case [11] were carried out, respectively.However, these studies mainly paid attention to the electromagnetic force itself, while the consequent responses such as the mechanical vibration deformations and strains, which will cause dynamic wearing to the winding insulation, have been rarely investigated in depth.What is more, in order to get a convenient analysis, scholars usually calculated the static integral force when the phase current is largest [12] as the most serious condition, whereas the moment when the phase current is largest may not be the exactly most serious time due to the negligence in the harmonic characteristics, which are sensitive to the mechanical vibration response.
Actually, it is more significant to find out the exact positions where the winding insulation endures the most serious mechanical response [13], including the static stresses, the strains, and the deformations (reflecting the vibration amplitude) [14].This is the basis for the winding design and manufacturing improvement [15].Reference [16] deduced the formula of the electromagnetic force density and the coupling vibration equations of the nonlinear magnetism-solid, based on which the amplitude of vibration by a numerical example is obtained.However, the serious vibration position of each end winding is not studied.
In view of the difficulty of 3D numerical solution due to memory limit, many researchers manually loaded a relatively stable electromagnetic force value to the 3D FEA model to solve the mechanical response on the end windings.Although this method can generally get the force distribution tendency, the result is not precise enough.The exact time for each winding to get the most critical stress and wearing threat requires transient electromagnetic-mechanical coupling response analysis.
As an improvement, in this paper, we propose a simplified model and carry out a fully comprehensive transient analysis on the electromagnetic force behaviour.The analysis includes not only the distribution regularities of the electromagnetic force, but also the force harmonic characteristics, and the mechanical response properties including the von-misses stress, the total deformation (reflecting the vibration amplitude), and the von-misses strains.

Study Object.
In this paper, we take the QFSN-600-2YHG type turbo-generator as the study object.The stator has 42 slots and the axial length is 6300-mm.The key parameters of the generator are shown in Table 1.
In the generator, double-layer lap windings are used, as shown in Figure 1(a).The three-phase windings are connected in double-wye forms.Since the winding distribution of the three phases is symmetric, the winding connection of Phase A indicated in Figure 1(b) is used as a typical illustration.The line part of the stator winding is laid in the slot, while the end part extends outside the stator core and forms a cone shape in two layers.The end winding is very complex and includes three parts, that is, the line part outside the slots, the involute part, and the nose part [1], as indicated in Figure 1(c).
To clarify each coil, the winding is indicated by the slot number of the line part on the upper layer.For example, Winding 1 indicates that the coil is composed of the upper layer line part in Slot 1, the bottom-layer line part in Slot 18, and the end parts connecting these two line parts.

Electromagnetic Force Analysis Model.
Since the line part of the winding is embedded inside the slots and bounded by closing wedges, the equivalent mathematical model of this part is a multi-supported beam.The main insulation wearing factor for this part is the static stress caused by the electromagnetic force and transferred from the stator core.Generally, this part is much less possible to be damaged because of its multi-restraining by the wedges and slots.In this paper, we primarily focus on the end part of the windings.
Since it is far more complex to calculate the magnetic field in the end zone precisely, the image theory is employed to transform the two-media problem into a one-media problem.Then, the image current obtained by the method is further used to deduce the magnetomotive force (MMF) of the stator windings.Finally, the electromagnetic force density is obtained based on the Ampere force formula.
Taking Phase A windings as an example, their image current is shown in Figure 2(a), and the generated electromagnetic field can be approximated as in Figure 2(b).Comparing Figure 2(b) with Figure 1(b), it is obvious that the approximate image current indicated in Figure 2(b) is the same as the current in the line part of the stator windings.Correspondingly, according to the relation between the rotor MMF and the stator MMF indicated in Figure 3 (in the figure   denotes the rotor MMF, while   denotes the stator MMF), the MMF in the stator end zone can be written as where   is the fundamental MMF produced by the image currents,  is the angle to indicate the circumferential position of the air-gap between the rotor and the stator, and Ψ is the internal power angle of the generator.Further, the endzone magnetic field can be assumed as a revolving field at the frequency of . = 2, and  is the mechanical rotating frequency of the rotor.A X The line part The nose part The involute part    Defining the circumferential position at the upper-layer line part of winding I as   , the induced current of this winding can be approximately described by where   refers to the angle between two slots, and 8.5  indicates the central position of the two line bars of the same coil (the pitch distance is of 17 slots).Λ 0 is the air-gap permeance per unit area,   is the MMF produced by the field windings,   is the exciting current of the generator,  is the turn number of the exciting winding for each pole,  is the axial length of the stator core, V =  ( is the inner radius of the stator core) is the rotating velocity of the magnetic field, and  is the reactance of the winding.
For the sake of clarity, we use a bar to indicate a space vector and a dot to indicate a time vector.Then the transient electromagnetic force density at an arbitrary Point  whose circumferential position is (  +   ) can be written as Feeding ( 1) and ( 2) into (3), the force density can be further expressed as where   is the angle between the current vector and the magnetic field density (MFD) vector at point .Then the general electromagnetic force on the whole bar can be obtained via the integral performance.
where  is the angle between the force density vector and the XY plane (the section plane of the stator core), and  is the circumferential angle that has the same meaning as previously mentioned.
As indicated in (4), it is evident that the frequency of the electromagnetic force is 2.Since the normal MMF contains each odd harmonics (the results shown in (4) only consider the 1st harmonic, which is the primary one), there should be other even harmonic components like 4, 6, 8 as well if considering the higher harmonics.Given that the amplitude will be deduced as the harmonic order goes higher, the electromagnetic force at 2 will be the major component and can be used as a characteristic component for the further study.

Mechanical Vibration Response Model.
The mechanical vibration of the stator winding can be treated as the periodic movement of the mass points (the periodic deformation of the winding).Mathematically, the periodic movement can be expressed as a displacement array of the mass points, whose kinetic equation can be written as where  is the mass matrix of the winding,  is the damping matrix,  is the stiffness matrix,   () is the exciting load matrix, that is, the electromagnetic force matrix, and () is the displacement response matrix of the structural mass points.Here the displacement response matrix is periodic and can be treated as the mathematical description of the mechanical vibration.
Generally, the physical model for the line part of the winding is a multi-supported beam because this part is all embedded in the stator slot and constrained by slot wedges and sealing, which have somewhat damping effects.However, for the end part of the winding, the physical model is an overhanging beam because this part is out of the slot.Consequently, mechanical vibrations will be primarily caused to the end part, and this part will be more probable to get an insulation wearing.

FEA on Electromagnetic Force
3.1.Models and Settings.Since the turbo-generator is very large, the meshing and calculating operation consumes too much memory that the computer spends too long a time to get a result.We use a simplified 3D FEA model, which only employs the end windings and one-tenth of the stator core (with the axial length of 630mm, while the actual length of the stator core is 6300mm).To verify the effectiveness of the simplified model, the 2D FEA calculation is carried out as well for a comparison.
The physical 2D and 3D FEA models for the QFSN-600-2YHG type turbo-generator are indicated in Figures 4(a stator windings and the exciting windings for both the 2D and 3D physical models are shown in Figures 4(e) and 4(f), respectively.
To obtain a more precise result, varied element length settings are adopted for different regions.The element length of the windings and the air-gap is set to the finest meshing grade, while the other parts are meshed with a coarser element length.The mesh results of the 2D and the 3D models are shown in Figures 3(b) and 3(d), respectively.There are 8710 elements for the 2D model and 197352 elements for the simplified 3D model, and the detailed mesh data is shown in Table 2.During the calculation, the vector potential boundary is employed, setting the outer stator edge (2D) and the outer stator surface (3D) as the potential reference.The rotating speed of the motion band is set to 3000rpm, and the exciting current is set to 4128 A, which is the rated value for the actual performance.

Analysis on Electromagnetic Force on Line Part.
Since the 2D model is actually a sectional model of the stator-rotor core and the line part of the windings, it cannot obtain magnetic force results for the end windings.Therefore, the further comparison between the 2D and the 3D models generally focuses on the results of the line part.
The integral forces on the line part obtained from the 2D calculation are indicated in Figure 5(a), while the result proposed in [17] is indicated in Figure 5(b).It shows that the force varying tendency of the 21 upper bars in this paper and in the reference paper is pretty similar.The difference in the exact force values is because the load of the generator in this paper is different from that in the reference paper (we cannot find the exact load data used in the reference paper).Generally, it is suggested from Figure 5 that the result obtained in this paper by the 2D calculation is valid.This offers a credible basis for the further comparison between the 2D model and the proposed simplified 3D model.
Since the current is a key factor to calculate the electromagnetic force acting on the winding (see ( 3)), we firstly check the current results obtained from the simplified 3D model.The currents of Phase A obtained from the 2D and the 3D models are indicated in Figure 6.they show that the  2D current is larger than the 3D result.The reason for this is that the 2D analysis does not consider the leakage magnetic field that will lose part of the conversion energy.Instead, all of the energy is taken into the current conversion during 2D calculation.However, it is still evident that the currents of Phase A obtained from the 2D and the 3D calculations match well with each other, which reflects the availability of the simplified 3D model in this paper.Given that the line part of the stator winding contains two segments, that is, the segment inside the slot and the segment outside the slot (see Figure 1(d)), it is necessary to further study the detailed electromagnetic force at different positions of the line part.In this paper, we set 4 probe points on the line part; see Figure 1(d).These four points have the same values in both X coordinate and Y coordinate, while the Z coordinate values of points 1-4 are 0mm, 40mm, 160mm, and 280mm, respectively.Since Point 1 is at the very position of the stator core edge in the axial direction, the electromagnetic force result of Point 1 is obtained from the 2D model, while the results of the other three points are obtained from the 3D model.
Typically, in a stable running condition, the air-gap magnetic flux density as well as the three-phase windings is symmetrically distributed during the rotor's rotation at a generally constant speed.Consequently, the electromagnetic force on the stator windings will be also symmetrically distributed in both time and space.In order to show the detailed force characteristics, winding 1 is taken as an example for the illustration.The electromagnetic force density variations of points 1-4 on the upper-layer line part of winding 1 (see Figure 1(d)) are plotted in Figure 7.
As indicated in Figure 7, the results from the 2D and the 3D models have the same time-variation law of sine function at the frequency of 2, which is consistent with (3).The electromagnetic force density amplitudes of points 1, 2, 3, and 4 are 3381, 3035, 3167, and 2178 KN/m 3 , respectively.The amplitudes of points 2, 3, and 4 are less than point 1 by 10.23%, 6.33%, and 35.58%, respectively.This is because points 2, 3, and 4 have a smaller and smaller distance to the end part where the magnet leakage effect exists (the MFD will be weaker in the end region).This phenomenon is also in accordance with the current result indicated in Figure 6.

Analysis on Electromagnetic
Force on End Part.Since the three-phase windings are symmetrically distributed, the electromagnetic forces on the three-phase windings should be similar as well.Limited by the space, only Phase A windings are analysed.
The electromagnetic force obtained from the 3D model on the end part of winding 3 when the phase current gets its max value (0.0854s) is indicated in Figure 8(a).As shown in the figure, the force vector has the largest amplitudes at two points, that is, the joint between the involute and the line part, and the nose-top.Along the axial direction, as it shifts away from the joint, the force amplitude will also be varied.
To get more detailed information about the electromagnetic force on the end part of the winding, the force amplitude varying tendency is indicated in Figure 8(b), where the data is quite similar with that in [2].In the figure, the original point is at the joint point.Since there are 7 coils in serial for one phase branch, it is necessary to obtain the max force amplitude for each coil.The maximum force density of Phase A windings is listed in Table 3, where it shows that Winding 7 has the largest force density, while Winding 3 has the smallest.
In addition to the force density, the integral force is also significant to the winding insulation wearing.The integral forces in time domain are indicated in Figure 9.It is shown that the forces are of sine-like function whose angular frequency is 2.The forces reach their peaks near the time t=85.4mswhen the current of Phase A gets its maximum    value.At this moment, the MMF produced by the image current rotates to the middle position of Phase A windings.Therefore, the phase current and the electromagnetic forces will get their max values.However, the moment for the force on each winding to arrive the peak is different, and it does not accord with the maximum current time exactly.Consequently, we cannot get the most serious states for each winding only by static analysis based on the maximumcurrent time.
To indicate the general force distribution on all of the coils, the integral force , and its components  x ,  y and  z respectively in X, Y, and Z direction, are shown in Figure 10 at the moment t=85.4ms.
As indicated in Figures 10(a)-10(d), the result illustrates a sine-function distribution tendency of the integral forces in the circumferential direction.In addition, there are two sine cycles for the resultant force and the Z direction component force, indicating that the frequency is 2.This result is consistent with (3).That means the force curves for windings 1-21 are the same as those of windings 22-42.This rule is the result from the opposite current as well as the magnetic field directions of the two-group windings.However, for the X direction and the Y direction forces, there is only one general sine-wave cycle. and  can both be decomposed in the radial and the tangential directions according to the following approximate equations: According to ( 5) and ( 7),  and  should be a square of cosine or sine.The square action will double the frequency of the curve, that is, cos 2 () = 0.5[cos (2) + 1].This explains why the resultant force has two general cosine-wave cycles, while the X direction and the Y direction force curves have only one general cosine-wave cycle.Meanwhile, it also indicates that the radial and the tangential force curves will

Analysis of the 2nd Harmonic Force on End Winding.
As previously analysed in Section 2.2, the electromagnetic force components are even harmonics, with the 2nd harmonic as the most prominent; see (4).In this paper, we also calculate the harmonic results of the electromagnetic force, and the force frequency characteristics of Winding 1 are shown in Figure 11.
As indicated in Figure 11, it is obvious that the 0-Hz component and the 2nd harmonic are the most conspicuous.The 0-Hz force component will only generate static stress on the windings, while the 2nd harmonic component will produce dynamic mechanical vibrations at 2.Comparatively, the static stress has much less threat to the insulation wearing than the dynamic vibration.In addition to the 2nd harmonic component, there are also other even harmonics like the 4th, the 6th, and so on.These pulsating components will also cause vibrations at their same frequencies.However, as the harmonic order goes higher, the amplitude of the force component will be decreased.Consequently, the vibration amplitudes caused by these components will also be decreased.That is the very reason why, in the practical monitoring, the 2nd harmonic of the force and the vibration is employed as the characteristic component.
To carry out a further study on the characteristic force component, the amplitudes of the 2nd harmonic of the integral forces for windings 1-42 are indicated in Figure 12, Mathematical Problems in Engineering where   denotes the resultant 2nd harmonic force, while    ,    and    refer to the components of the 2nd harmonic force, respectively, in X, Y, and Z direction.
It is shown in Figure 12 that there are six sine-like cycles for the 2nd force.The reason for this is that there are three phases, and each phase has two parallel branches; see Figure 4(e).Moreover, the 2nd harmonic forces on windings 1-21 are the same as those on windings 22-42.This is because these two groups of windings, which are actually the two parallel branch winding groups, are always on the opposite side in the symmetrical space distribution and connected in the opposite direction; see Figure 1(b).When the N-pole of the magnetic field is acting on one group of windings, the S-Pole will be acting on the other group so that the current and the force of these two winding groups will be the same.

Study on Mechanical Responses
Since the electromagnetic force will cause significant mechanical response including vibrations, deformations, and equivalent stresses and strains, the winding insulation is probably damaged after standing these mechanical responses for a long time.The key point to protect the winding insulation is to obtain the dangerous positions where the winding endures the most serious total deformations, equivalent stresses, and elastic strains, so that proper measures can be carried out.
As previously analysed in Section 2.3, the insulation of the line part will be much safer than the end part due to the structure difference.Therefore, in this paper, we mainly calculate the mechanical responses of the end winding.The winding material is defined as copper alloy whose Young's Modulus is 1.1×10 5 Mpa, Poisson's Ratio is 0.34, and yield strength is 280Mpa.
During the calculation, the electromagnetic force density data obtained from the 3D model in the electromagnetic module is imported into the structural module to further calculate the mechanical responses, as indicated in Figure 13.The line part is constrained by a fixed support.Automatic hexahedral solid element is adopted and 109745 nodes are generated for each end winding.The maximum and the minimum equivalent stresses, the total deformations, and    the elastic strains of the 7 windings, when Phase A current is on the top, are indicated in Figure 14, while the detailed largest responses of each winding are illustrated in Figure 15.
As indicated in Figure 15, the max mechanical response values as well as the moment to get the max values will shift from one winding to another.Generally, both the max stress and the max elastic strain appear on winding 1 at 0.015s, while  the max total deformation appears on winding 1 at 0.0195s.The max stress is 131.97MPa, which does not exceed the yield strength limit of the winding material.When the current is largest, the max stress is 129.5MPa.This means that the very time when the current gets the peak value is not the exact moment for the winding to endure the max stress.Therefore, the max response cannot be obtained by the easy analysis based on the largest phase current moment.According to Figure 15, the max deformation is 1.39 mm.The practical data may be even larger, because the analysed winding model is established as ideal solid, and the rigidity is larger than the actual winding.
As indicated in Figure 14, comparatively, both the nosetop and the middle of the involute have the most intensive stresses, deformations, and elastic strains.In addition, it is also suggested from Figure 15 that the line-involute joint has the largest stress and strain, while the nose-top and the middle of the involute have the largest deformation.
On one hand, the severe stress will intensify the degradation of the insulation material.On the other hand, the intensive vibration amplitude, which is reflected by the deformation, will increase the friction of the insulation material and will make windings impact against each other.Subsequently, after a long period, the fatigue effect will act, and crevices and holes will probably be caused at the nosetop, in the middle of the involute, and at the line-involute joint.
Practically, we have also found actual damage at these mentioned three positions, especially at the nose-top.Figure 16 illustrates the damage photos of different generators, and the damaged districts follow well the dangerous locations concluded in this paper.Correspondingly, these three parts need special reinforced measures during the manufacturing, the assembling, and the monitoring process.
To illustrate a more visualized result about the mechanical response distribution principles, we calculate the minimum and the maximum values of the transient stresses for windings 1-7, as indicated in Figure 17.It is shown that Windings 1 and 7 have both the least minimum stress and the largest maximum stress.That means, the windings between two phases will endure more serious mechanical responses than others and therefore need more attention during both design and monitoring periods.The calculation data of Phase B and Phase C also confirms this phenomenon.
In addition to previous findings, there is still another interesting phenomenon, which can be obtained by comparing Figure 9 with Figure 15, and Table 3 with Table 4.It is shown that the winding, which endures the maximum mechanical response, is neither always consistent with the one that has the largest resultant electromagnetic force nor directly the same as the winding that affords the most intensive electromagnetic force density.This is significant because it is different from people's assumed recognition.The very reason for this phenomenon is that the force density is unevenly and nonlinearly distributed on the involute.Taking Figure 18 as an example, case 1 and case 2 have the same integral force 7N, while case 1 has even larger force at the point P.However, case 2 has the larger bending moment to the point P than case 1. Subsequently, the equivalent stress and deformation on point P in case 2 is larger than that in case 1.

Conclusions
In this paper, a comprehensive analysis on the transient electromagnetic force behavior of stator windings in a QFSN-600-2YHG type turbo-generator is presented.The work is generally based on a proposed simplified model and the 3D finite element coupling calculation.The simplified model uses the approximate solution of the image current to conveniently obtain the analytic formula of the electromagnetic force on the end windings, and employs only the end windings   and one-tenth of the stator core for the 3D FEA to save the calculating memory, but meanwhile obtain the qualified electromagnetic forces as well as the mechanical responses.The result comparison between the proposed model and the 2D model as well as the data in published references confirms the validation of the simplified model.This will be beneficial to improve the transient calculation speed in large turbogenerators.
It is shown that the electromagnetic force has only even harmonics, with the 2nd harmonic being the prominent component.Under the action of the electromagnetic force, the nose-top, the joint between the line part and the end part, and the middle of the involute are the three most dangerous positions for the insulation wearing.This finding can be used as a reference to improve the insulation manufacturing and assembling techniques as well as the winding monitoring performance.
What is more interesting and significant, it is found that the winding, which endures the maximum mechanical response, is neither always consistent with the one that has the largest resultant electromagnetic force nor directly the same as the winding that affords the most intensive electromagnetic force density.This is different from people's assumed recognition and therefore needs special caution during the design and the maintenance processes.

Figure 1 :
Figure 1: Structure of stator windings: (a) sectional view of three-phase winding distribution, (b) winding connection of Phase A, (c) structure components of the end winding, and (d) probe point set for further analysis.

Figure 5 :Figure 6 :
Figure 5: Comparison of integral force on line part: (a) 2D result in this paper, and (b) result proposed by Ref. [17].

Figure 7 :
Figure 7: Electromagnetic force density on the line part.

Figure 8 :
Figure 8: Electromagnetic force density on the end part of winding 3: (a) linear density from reference; (b) volume density from 3D FEA.

Figure 9 :
Figure 9: Integral force of Phase A windings.

Figure 11 :
Figure 11: Harmonics of electromagnetic force on winding 1: (a) resultant force, (b) component in X direction, (c) component in Y direction, and (d) component in Z direction.

Figure 12 :
Figure 12: 2nd-harmonic amplitude distribution: (a) resultant force, (b) component in X direction, (c) component in Y direction, and (d) component in Z direction.

Figure 13 :
Figure 13: Support and imported load settings: (a) support, (b) imported body force density.

Figure 16 :Figure 17 :Figure 18 :
Figure 16: Real insulation damage photos in different generators: (a) nose-top damaged in a 200MW generator (the black area indicates the insulation burned districts), (b) nose-top damaged in a 15MW generator, (c) the middle of the involute part damaged in a 20MW generator, and (d) joint part damaged in a 300MW generator (caused by overstress-induced spacer deformation and friction).

Table 1 :
Key parameters of the study object.

Table 3 :
Maximum force density of phase A windings.

Table 4 :
Maximum stresses and deformation statistics on Phase A windings.