A Model for Stator Eddy Current Losses Due to Axial Flux in Synchronous Generators at Steady State and under Load Angle Oscillations

Abstract The load angle of a synchronous generator connected to a power grid has an eigenfrequency that depends on the operating conditions. The existence of an eigenfrequency can make the generator sensitive to electrical and mechanical disturbances and motivates the use of damper windings and power stabilizing systems. Another reason for keeping load angle oscillations small is that they increase the eddy current losses in the stator core and clamping structure due to axial leakage flux. This is shown with transient finite element analyses and can be explained by a loss model derived from a phasor model of the eddy current loss density.


I. INTRODUCTION
For each point of operation, a synchronous generator connected to a power grid has an eigenfrequency with which the rotor oscillates relative to the stator current contribution of the magnetic field. Mechanical or electrical noise with a wide spectrum and a low intensity can trigger and maintain the oscillations. The natural frequency is often in the range 1-3 Hz [1,2]. Transient stability of the synchronous machines and demands on keeping the voltage and frequency on the power grid within tight limits motivate the use of rotors with high inertia, damper windings and Automatic Voltage Regulation (AVR) with Power System Stabilizer (PSS). Control systems for improved transient stability at large disturbances have been suggested [3]. Also High Voltage Direct Current (HVDC) has been suggested for transient stability [4,5]. For generators connected to a strong grid, the incentive for keeping load angle oscillations small can increase if it can be shown that load angle oscillations lead to increased losses. The importance of axial leakage flux for end region heating was recognized in 1920 [6] and has been in focus many times. The effect of steady state operating conditions on axial magnetic leakage flux density and losses in the stator steel of a synchronous generator has been investigated [7,8]. The influence of pole slipping on axial magnetic leakage flux has been measured [9]. That eddy current losses due to axial leakage flux increase at load angle oscillations has not been shown previously. Furthermore, when eddy current power losses due to axial leakage flux are high, other iron power losses can also be expected to be high, at least at no load [8].
The purposes of this paper are to derive a loss model and to show the effect of load angle oscillations on losses caused by eddy currents in the stator core of a synchronous generator. Because of limited computer resources, Finite Element Analyses (FEA) of oscillations of a machine with small poles have been done. Since the eddy current power losses in this machine were found to be low both absolutely and relative to the rated power, simple supplemental steady state FEA have been done for estimation of the effect of the machine size and the stator core permeability. Furthermore, since the damper windings have not been included in the FEA with loss calculations, a system of ordinary differential equations has been solved to calculate the effect of damper windings on winding currents.
ANSYS Design Modeler and time stepped FEA with ANSYS Maxwell [10] have been used for machine design and calculation of inductances and losses in this work. All calculations not done with FEA have been done with the programming and calculation software MATLAB [11].
In this paper, a bar above a letter denotes a complex quantity. A bold symbol denotes a vector or a matrix.

Phasors, Angles, Voltage Equations and Equations of Motion
Since the natural frequency of a large synchronous generator is much lower than the electrical frequency of the power grid, the phase voltages and phase currents vary approximately harmonically but with time dependent amplitudes during load angle oscillations. The almost harmonic variations make phasors appropriate for description of the electrical state of the generator. Phasors and angles used in this article are shown in Figure 1 with an infinite bus rotor reference system with direct axis d IB and quadrature axis q IB . The root mean square (rms) terminal phase voltage is V , the rms stator current is I , the field current is i f ¼ Ài f , the rms electromotive force (emf) no load phase voltage is E ¼ jE, and the rms phase voltage at the infinite bus is V IB : There is normally a transformer between the generator and the infinite bus. In that case, V IB is the infinite bus voltage referred to the generator side of the transformer. It is assumed that the transmission lines are so short between the generator and the infinite bus that the current is approximately the same along the line, i.e. I IB ¼ I : Prior to any disturbance, the generator is assumed to be synchronized with the infinite bus reference frame. Subscript 0 denotes the steady state starting position for load angle oscillations. The electrical rotor angle is h from the stator phase a axis to the q axis, the time is t, the load angle is d ¼ argð EÞ À argð V Þ, the load angle from the infinite bus phase voltage phasor to the emf phasor is d IB ¼ argð EÞ À argð V IB Þ, and the phase angle of the stator current is a from the d axis. The change in d IB is b ¼ d IB À d IB, 0 : The rotor electrical angular speed is x ¼ x s þ db=dt where x s is the synchronous and infinite bus angular speed.
The steady state stator voltage equation is where X d is the d axis stator reactance, X q is the q axis stator reactance and R is the stator resistance per phase, I d ¼ I d ¼ I cos a and I q ¼ jI q ¼ jI sin a: During load angle oscillations, (1) is approximately fulfilled. Normally, the variation of pole flux linkage is smaller than the variation of i f during rotor angle oscillations [2]. In such cases, a more convenient voltage equation is the constant flux linkage model [2,12], where E 0 is the q axis emf which is proportional to the pole flux linkage, and X 0 d axis reactance where L df ¼ -ffiffiffiffiffiffiffi ffi 3=2 p M f is the mutual inductance between the field winding and the fictitious stator d axis winding, M f is the maximum of the mutual inductance between a phase winding and the field winding, L f is the self inductance of the field winding, X 0 td ¼ X 0 d þ X tr , X tq ¼ X q þ X tr , and R t ¼ R þ R tr : Here, X tr is the reactance and R tr is the resistance of the combination of a transformer and a transmission line between the generator and the infinite bus.
Most general and accurate is to use transient voltage equations with damper windings like in [13]. For each pole in a generator called G2, there are four electrically insulated, unbranched damper windings distributed such that each pole is symmetrical with respect to the d-z plane. The symmetry makes it correct to include only one pole in a FEA. Generator details are given below in Section 2.4. Unbranched windings have been used since their inductances can relatively easily be calculated as functions of the rotor angle in FEA. The inductances have been used in solution of voltage equations and equations of motion with solvers for ordinary differential equation systems. The transient voltage equations for the phase windings a, b, c, the field winding f and the four damper windings, U, D, W and Q are q for P to be orthogonal. In this case, with inductance matrices L S , M and L R , and exponent T denoting a transpose, the transient voltage equations formulated as equations for the current derivatives in the rotor reference frame can be written as d dt where 0 is a square matrix of zeros, The dq0 transformation with k ¼ ffiffi 3 2 q and balanced phases The swing equation of motion of the rotor of a synchronous generator is given by [14] J dx m dt ¼ T m À T e (12) where J is the polar mass moment of inertia, T m is the mechanical torque on the rotor, T e is the electromagnetic torque on the rotor, x m ¼ x=p is the mechanical angular speed where p is the number of rotor pole pairs. If the damper currents are not calculated and not considered in T e , their damping can be included in a separate electromagnetic damping term. In this special case, (12) can be written as where x m, s ¼ x s =p, and D 0 is a damping coefficient. The load angle oscillation analyses are restricted to small variations of rotor angular speed such that x=x s % 1: Therefore, multiplication of (13) by x m gives where D ¼ D 0 =p, P m is the mechanical power on the rotor shaft, P e is the electric power developed by the generator. The electric power can be written as where P IB is the power delivered from the generator to the infinite bus. In case of steady or almost steady operation, by use of (2), P e can be expressed as For numerical solution, (14) can be written as a system of two equations that, with use of In the case of small b and db=dt, (17) can be linearized around b ¼ 0. That way an equation of motion of the same form as for a mass point connected to a damper and a spring with constant stiffness can be obtained [2,15]. The period of load angle oscillation, t p , can be estimated from an expression for the period of a torsional pendulum, t p ¼ 2p ffiffiffiffiffiffiffiffiffiffi ffi J=K p , where K is the torsional stiffness [15]. For a generator connected directly or indirectly to an infinite bus, the period can be estimated by where ðdT b =dd IB, m Þ av is the average torsional stiffness, the so called synchronizing coefficient [1], T b ¼ ðP e À P m Þ=x m % ðP e À P m Þ=x m, s is the electrical torque minus the mechanical torque on the rotor, and d IB, m ¼ d IB =p: Furthermore, P m is independent of d IB . The period can therefore be written as where H is the inertia constant and S is the nominal apparent power of the generator. With damping and resistance neglected, (15), (16) and (19) give If s is the length scale and if a generator is scaled with the proportions and the peak flux density maintained, J is proportional to s 5 , flux and voltage are proportional to s 2 if resistance is negligible, and current and inductance are proportional to s. Then the period according to (20) increases in proportion to s.

Phasor Models of Magnetic Field and Eddy Current Density
The starting point is the assumption that the magnetic flux density, B, and magnetic field strength, H, change harmonically with time. A relation assumed to hold in this study is where the permeability l i is constant, and i ¼ x, y or z. A simple phasor model is where i ¼ x, y or z. The model (22) is a generalization of the phasor model for the axial component, B z , in [7]. The parameters k f , i , k d, i , k q, i , f , i , d, i and q, i depend on frequency, permittivity, permeability, conductivity, geometrical parameters and the location and are assumed to be independent of i f , I d and I q : The coefficients k f , i , k d, i and k q, i include damping and can be chosen to be positive because of the exponential factors. With r denoting conductivity and e denoting permittivity, the condition r ) xe is well fulfilled for steel at the fundamental frequency in a power generator. Then Amp ere's law is simplified to where J is the current density. Equation (22) inserted in (23) gives that the components of J can be written in the same form as (22), i.e.
The square of the magnitude of (25) is

A Model for Eddy Current Loss
The time average of the eddy current loss density during a cycle with period 2p=x at steady state and approximately at load angle oscillations can be written as Volume integration of (27) with use of the mean value theorem [16] for a part with volume V and constant conductivity gives that the total time averaged eddy current loss during a cycle is where subscript me refers to the mean value in the volume. In a linear generator and approximately in a rotating machine with many poles, the fundamental frequency component of the eddy current density at steady state can be approximated by J ðx, y, z, tÞ ¼Ĵ ðy, zÞ cos px where x, y and z are Cartesian coordinates of which x is in the direction of motion of the translator or rotor pole, / is a phase angle, and s p is the pole pitch. In the special case when V spans a pole pitch and therefore a period of J ðx, y, z, tÞ 2 , Equation (29) gives that the spatial average in the x direction of the loss density is equal to the time averaged loss density, i.e. p ed, me ðy, zÞ ¼ 1 The instantaneous eddy current loss at steady state in a stator part with volume V spanning a pole pitch is therefore the time averaged loss, i.e.
In case of load angle oscillations, also the amplitudeĴ and P ed, av are time dependent but so slowly that (31) is approximately valid. The time average of P ed, av in a b period t p is The losses calculated at different points of operation with (31) can be used for estimation of the coefficients in (28). Equation (28) can be simplified to Each point of operation gives a set of values of P ed, av , a, b, c, d, e and f that can be used in (33) as a condition on the unknown parameters x, y, z, u, v, and w. These are components of the solution vector of the equation system

Finite Element Models, Parameters and Materials of the Generators
In this work, short salient pole synchronous generators, G1 without and G2 with damper windings, have been designed and modeled in ANSYS Design Modeler only for the study of load angle oscillations. The analyzed machines exist only as finite element models designed to be so simple and small that they can be used for time stepped oscillation analyses in the available computer. For symmetry reason, it is sufficient to include in FEA one half of the axial length of one pole pitch of the studied generators. Figure 2 shows G1 and G2. Figure 3 shows the mesh of the stator core and clamping structure of G1. The mesh is good enough to capture the trends predicted by the loss model in (33). The rms element side lengths are 4.5 mm in the core, 3.9 mm in the fingers and 4.9 mm in the ring within a radius of r lim ¼ 1565 mm. Each stator tooth and clamping finger has been cut off at an angle of 45 from a plane perpendicular to the z axis along the shaft. Geometrical parameters, inductances and resistances of the simulated generators are given in Table 1. Material data and the numbers of finite elements are given in Table 2. Diagonal permeability and conductivity matrices have been used. Only the rotor rim and poles with damper windings disregarded have been included in the calculation of J: Figure 4 shows unbranched, electrically insulated damper windings of G2. After preliminary analyses with electrical contact between the clamping ring and the clamping fingers, insulating boundary conditions have been used on the clamping fingers and between damper windings. For resistance calculation, the field, damper and stator phase windings have been assumed to be made of annealed copper with conductivity 48.6 MS/m at a temperature of 70 C.
Generators G3 and G4 have been used only for the analyses of the effects of the machine size on losses at steady state. Generator G3 differs from G1 by having a longer and laminated rotor core and partly other materials. Generator G4 is G3 upscaled a factor 1.2 in every direction. Within the radius r lim ¼ 1565 mm in G3 and r lim ¼ 1878 mm in G4 in the stator core, the relative permeability is 100 in the x and y directions and 9.3262 in the z direction, corresponding to a stacking factor of 0.9. The lower permeability in the teeth is realistic at flux densities of around 1.9 T according to extrapolation of data for nonoriented steel M350-50A [17]. The rotor core is axially equally long as the stator core. In addition, eddy currents in G3 and G4 have only been calculated in the clamping structure and the end of the stator core so that a slightly finer mesh could be used. The core end length with eddy currents has been chosen to be 8% of the core length, i.e. 12 mm for G3 and 14.4 mm for G4. Finite element size specifications are the same for G3 and G4 in the air gap and the eddy current carrying parts except the outer part of the clamping ring. The rms element side lengths are 3.2mm in the core end, 3.5 mm in the fingers and 4.9 mm in the ring within a radius of 1565 mm in G3 and within a radius of 1878 mm in G4.
In FE models spanning only one pole pitch, antiperiodic boundary conditions have been used on boundary planes at a constant cylindrical angle coordinate. At the other boundaries, the normal component of the magnetic field strength is zero by default.  In machine guidelines [18], the amplitude 0.85-1.05 T in the air gap is recommended for the fundamental frequency in a salient pole synchronous machine. The spatial distribution of the radial flux density, B r , can be written as a Fourier series [16] B r ð#Þ ¼ where # is an electrical angular position along the air gap in the direction of rotation, and the Fourier coefficient is The amplitude of the fundamental component of B r in the air gap isB gap, 1 ¼ 2jc 1 j: However, an estimation,B      Table 3. The generators have been assumed to be directly connected to an infinite bus. After an iterative design process aiming at a machine design that givesB 0 gap, 1 % 1 T at steady state with power factor 1 and reasonable currents, the steady state I was chosen to be 2000 A for G1. With the phase voltage chosen to be V ¼ V IB ¼ 81 V, and cos u ¼ 1, (1) gives E % 101:05 V. The same E was calculated with no load analysis in Maxwell 3 D with i f ¼ 220.44 A. With these V and I, and E calculated with (1), i f at overexcitation and underexcitation was calculated from i f ¼ Eð sin u ¼ 0.6) i f ð cos u ¼ 1)/Eð cos u ¼ 1) and i f ¼ Eð sin u ¼ À0.6) i f ð cos u ¼ 1)/Eð cos u ¼ 1), respectively. Also for G3, the choice of I was 2000 A. The procedure for obtaining suitable V and i f for G3 began with a no load analysis with inductance calculation but no eddy currents. The field current 180 A gave E ¼ 74.11 V in FEA. With the inductances d-q transformed, (1) was used for finding V that gives E at sin u ¼ À0.6. A subsequent preliminary load analysis gaveB 0 gap, 1 ¼B 0 gap, 1, prel which was less than the desired 1.0 T. Therefore, V was scaled by the ratio 1.0 T/B 0 gap, 1, prel : The corresponding E ¼ 83.11 V was calculated by (1), and i f was calculated by 83.11/74.11Á180 A. For G4, the currents I and i f were taken from G3 and scaled up with the geometrical length scale factor s ¼ 1.2. A no load, FEA gave the corresponding E and the inductances. Finally, (1) was used for finding V that gives E. The value of V is somewhat higher than what would be obtained by just scaling the voltage of G3 by s 2 since the resistance is not negligible for the studied small machines.
The stator inductances in the

Load Cases
The initial states for the load angle oscillation analyses are given by Table 3.
With G1, the studied cases have either constant i f or constant v f . In cases of constant i f , (17) has been solved with constant P m and D ¼ 0. In each time step in the solution of (17), d IB ðbðtÞÞ is updated, V d is calculated by V d ¼ V sin d IB , V q by V q ¼ V cos d IB , I d and I q are calculated by (1), and P e by (15). A zero-initial condition on b and a  nonzero initial condition on db=dt that gives a desired b max > 0 have been used to simulate steady b oscillations with constant i f . Steady state has been simulated with initial With G1 and G2 in cases of constant v f , (17) has been solved with respect to b and db dt simultaneously with (7) with respect to the currents. In these calculations P m is where P m, 0 is the steady state initial value of P m , P is a pulse equal to 1 for t 1 < t < t 2 , 1/2 for t ¼ t 1 and t ¼ t 2 , and 0 otherwise, j is a scale factor for the disturbance amplitude, and t r is the resonance period found by repeated solving. The chosen values of j and D are high enough for b oscillations with approximately constant and desired amplitude, 12 , to occur within a 28 s duration of the disturbance. With initial power factor 0.8, constant v f and b max ¼ 12.0 , the underexcited case has been simulated with D ¼ 650 Nm, j ¼ 0.01434, and t r ¼ 190.4 ms for G1, and with D ¼ 0, j ¼ 0.025844, and t r ¼ 190.2 ms for G2. The overexcited case has been simulated with D ¼ 1200 Nm, j ¼ 0.01268, and t r ¼ 157.7 ms for G1, and with D ¼ 0, j ¼ 0.017001, and t r ¼ 157.5 ms for G2. For G2, D ¼ 0 since the damper windings are coarse enough to allow stable forced load angle oscillations at resonance frequency.
The effect of machine size on eddy current power losses due to axial leakage flux has been investigated with G3 and G4 at underexcited steady state with power factor 0.8. Figure 5 shows that the time averages of the eddy current losses in the whole stator core, all clamping fingers and both clamping rings at constant i f and P m increase with b max . Figure 5 also shows that the losses increase with decreasing excitation, which agrees with results of measured and FE simulated axial flux density [7,8]. The time averaged eddy current losses in the two studied cases with constant v f and no damper windings differ by less than one per mille from those at the same b max at constant i f . Figure 6 shows the phase angles as functions of time in the two cases with constant v f . The phase angles change similarly at constant i f . As expected from (20), the overexcited state has a shorter period than the underexcited. Figure 7 shows I d , I q and i f as functions of time in the cases b max ¼ 12 and initial power factors 0.8 underexcited and overexcited. The variation of I d is larger at constant v f than at constant i f , but since I d and i f increase simultaneously and give (obliquely) opposing magnetic field contributions when I d is positive, it is not surprising that the cases constant i f and constant v f give similar eddy current losses in the stator. With d defined by Figure 1, the stator current components change with d roughly according to (1). The peak in I d at d min in the overexcited state occurs also when i f is constant and is a consequence of a non-negligible R. A comparison between Figure  6 and Figure 7 shows that I q changes approximately in phase with d and b. So do P e and the apparent power.

Results with Generator G1
Twenty-one different points of operation with constant i f and db dt ¼ 0 have been used for numerical determination of the six model parameters in (33). The parameters are approximately x ¼ 35.8, y ¼ 3.69, z ¼ 4.19, u ¼ À15. 8, v ¼ À0.772, for the ring. With these parameter values, (33) has been used in the cases at constant v f for comparison between FEA and the loss model, as shown in Figures 8 and 9. Ripple and minor noise in the loss calculated in the FEA have been smoothed out. The bad core loss fit the first tens of ms is due to the fact that the eddy currents start from 0 in FEA, and it takes time for the eddy currents to build up. Figure 10 shows that the magnitude of the eddy current density is high in the end region. Preliminary FEA showed that the losses in the clamping structure increased if the electrical insulation on the clamping fingers was removed. Figure 11 shows the currents in the damper circuits. The damper winding currents are not only low but approximately 90 out of phase with I d , I q and i f in Figure 7. Comparisons between G1 and G2 show that the damper currents changed the maxima and minima of I, I d and I q by less than 4% except the maximum of I d at underexcitation that decreased by 16% from 211 A. The maxima and minima of i f changed by at most 0.2%.

Results with Generators G3 and G4
The stator eddy current power losses at underexcited steady state at power factor 0.8 andB

Results of Decreased Core Permeability at Steady State
At steady state with underexcitation, the change of relative permeability from 1000 to 100 in the x and y directions of the whole stator core in G1 resulted in more than doubled loss in the core alone and almost 7 times higher total average eddy current loss due to axial flux in the core and clamping structure of the studied stator although i f was not increased to compensate for the higher reluctance. Currents at load angle oscillation in four cases with load angle amplitude 12 . In the two cases with constant field voltage, the field current is shown scaled up for improved visibility. In the other cases, the field current is according to Table 3. Meanings of subscripts:

DISCUSSION
In generator G1 with permeability according to Table 2, the total average eddy current loss due to axial flux at steady state is only about 0.1 per mille of the rated power. However, this is very dependent on the length scale s. Whereas for a given flux density amplitude, the rated power is proportional to s 3 , the eddy current loss power is at least approximately proportional to s 5 according to Section 3.3. The latter is supported by the case where a harmonically time varying magnetic field, B(t) ¼ B 0 sin xt, in the axial direction of a solid conducting cylinder is approximately uniform over the circular cylinder cross section. In this case, the eddy current power loss is approximately where r is the radius, h is the axial height, and V is the volume of the cylinder. However, the scaling is complicated by both segmentation of the core and the clamping ring. In addition, the results in Section 3.4 show that a decreased permeability in the x and y directions of the stator core can increase the eddy current losses at least if the laminate efficiently limits the eddy currents caused by flux parallel to the laminate planes. The increase of the eddy current losses can then be explained by a deeper penetration of the axial flux into the core [19]. The loss model can be expected to work also for nonsalient pole generators. According to Figures 8 and 9, the loss model fits the loss data well except in the stator core at b min at overexcitation when the losses are at their minimum. The cause of this is yet unknown, but FEA shows that the radial component of the eddy current loss is relatively far from harmonically time varying in the stator tooth sides near b min at overexcitation.

CONCLUSIONS
At steady state and at load angle oscillations, it is possible to estimate the eddy current losses in the clamping structure and the stator core with a simple loss model of quadratic functions of the field current and the d-q components of the rms stator current. The relative sizes of current components and the parameters in the loss model (33) determine how important the loss terms are. The squares of the two stator current components are roughly equally important for the losses. The field current is more important for the eddy current losses in the clamping structure than in the stator core. This can at least partly be explained by the reluctances of the flux paths through air. The average flux path length through air at the generator ends is longer at high field current than at low field current whereas the average flux path length in the radial air gap is relatively independent of the point of operation according to flux density vector plots in FEA. Since the d component of the stator current is positive at overexcitation and negative at extreme underexcitation, the negative sign of the coefficient in front of the product of the field current and the d component of the stator current is in agreement with previous findings of higher axial magnetic leakage fields [7,8] and higher losses [19,20] at underexcitation than at overexcitation. The last two terms in the loss model (33) are relatively small. This is in accordance with the fact that the magnetic field from the q component of the stator current  Since the rms stator current increases with the load angle, the two terms containing the squares of the stator current components in the loss model can explain why the time averaged eddy current loss increases with the amplitude of the load angle oscillation.
For a given flux density in the air gap, the eddy current losses caused by axial leakage flux are approximately proportional to the length scale raised to 5 whereas the rated power is approximately proportional to the length scale raised to 3 if the resistance is negligible. Under such conditions, the importance of the axial leakage flux increases rapidly with the machine size.
he was employed at BAE Systems Bofors. There he worked with various technical calculations and writing technical reports. He also did technology monitoring in pulsed power and electromagnetic launch systems. Most of the assignments were within electromagnetics and heat transfer. In 2011-2017 he was a Ph. D. student at Uppsala University doing research related to axial leakage flux in synchronous generators. After receiving his Ph. D. degree in 2017, he 2018 wrote educational material on synchronous generators for hydropower. Since 2019 he is a consultant responsible for high voltage electrical requirements of climate systems on Volvo buses. His research interests include electric power components and systems.
Urban Lundin received the Ph.D. degree in condensed matter theory from Uppsala University, Uppsala, Sweden, in 2000. After a post-doc period he joined the Division for Electricity, Uppsala University, in 2004 where he is currently a Professor in electricity with a specialization towards Hydropower, where he heads the Hydropower group. He has been involved in the industrial implementation of research projects including the electrical spillway, magnetic bearing, and an excitation system to remove unbalanced magnetic pull. His research interests include synchronous generators, excitation systems and their interaction with mechanical components and the power system.