Coupled Electromagnetic–Fluid–Thermal Analysis of a Fully Air-Cooled Pumped Storage Generator Motor

: With the continuous increase in the capacity of the pumped storage generator motor, the overheating of the rotor area is becoming increasingly severe, which has a signiﬁcant effect on the safe and reliable operation of the machine. The heat dissipation of the machine rotor by fully air-cooled is one of the key technologies to develop the new generation of pumped storage generator motors. In this paper, the electromagnetic ﬁeld and ﬂuid–thermal coupled ﬁeld of a pumped storage generator motor are analyzed. The 2D transient time-step ﬁnite element model of the electromagnetic ﬁeld of a pumped storage generator motor is established, and the eddy current loss of damping bars of the rotor is calculated by the ﬁnite element method. The additional loss of the rotor pole surface is calculated by analytical method. The mathematical and geometric models of the 3D ﬂuid–thermal coupled ﬁeld of the pumped storage generator motor are established and calculated. The complex ﬂuid velocity distribution and the temperature distribution at different positions of the rotor under fully air-cooled fanless cooling conditions are investigated in detail. The calculated temperature of ﬁeld winding is compared with the measured value, and the result shows that the calculated result coincident well with the test data. This research provides the technical reference for the development and temperature rise calculation for large pumped storage generator motors.


Introduction
With the rapid growth of the global economy, the continuous excessive consumption of non-renewable energy has resulted in a series of environmental problems.The world is actively promoting the green energy transformation and promoting the highquality development of energy.Energy storage systems are vital electric grid facilities for promoting global clean energy transformation [1,2].Pumped storage is the largest-scale grid-connected energy storage technology.A pumped storage unit can perform tasks such as peak shaving, valley filling, phase modulation, and frequency modulation in the power grid, greatly promoting the integration of wind and solar power, and ensuring the safe and reliable operation of the power system [3].With the continuous development of pumped storage units, the power rating of pumped storage generator motors is increasing.The electromagnetic losses, field current, and temperature rise of the machine rotor increase obviously, which leads to overheating issues.Overheating of the machine rotor will seriously influence the safe and reliable operation of the pumped storage generator motor, which will cause instability of the power grid.
The fully air-cooled fanless cooling system of large hydrogenerators offers the benefits of robust reliability, low cost, and convenient operation and maintenance.However, the fully air-cooled fanless cooling system is an indirect cooling approach, and its cooling flow rule of cooling air and the temperature distribution of the machine rotor are analyzed in detail.

Establishment of the 2D Transient Electromagnetic Model
The 2D transient electromagnetic model established in this paper is based on a 145 MW (generating operation)-170 MW (pumping operation) pumped storage generator motor.The geometric dimensions and electromagnetic parameters of the studied machine in this paper are summarized in Table 1.The established 2D transient electromagnetic model is shown in Figure 1.The solving region of the built model includes a stator core, armature winding, air-gap, a pair of rotor poles, rotor yoke, field winding, and damping winding.The damping bars of one damping winding are numbered from 1 to 8 along the clockwise direction and electrically connected at the end.According to the machine's rotation direction, one side of the rotor pole is defined as the windward side, while the other is defined as the leeward side.A total of 21,980 triangle elements are used to divide the solving region and the mesh in the air-gap region and damping windings are refined to enhance the calculation precision.
To simplify the analysis, the following assumptions need to be proposed: 1.
The electromagnetic model is built as a 2D electromagnetic field and only the axial component of magnetic vector potential A and current density vector J are considered; 2.
The leakage flux outside the machine core is ignored.
In the solution region, the governing equations of the transient electromagnetic field and associated boundary conditions can be derived as follows: where σ is electric conductivity (in S/m), µ is permeability (in H/m), t is time (in s), and A z and J z are the magnetic vector potential and current density components along the axial direction.
Machines 2023, 11, x FOR PEER REVIEW 4 of 21 1.The electromagnetic model is built as a 2D electromagnetic field and only the axial component of magnetic vector potential A and current density vector J are considered; 2. The leakage flux outside the machine core is ignored.In the solution region, the governing equations of the transient electromagnetic field and associated boundary conditions can be derived as follows: where σ is electric conductivity (in S/m), μ is permeability (in H/m), t is time (in s), and Az and Jz are the magnetic vector potential and current density components along the axial direction.

Electromagnetic Field Calculation
The pumped storage generator motor has two working conditions: generating and pumping, which can be treated as a salient-pole synchronous generator and a salient-pole synchronous motor, respectively.Figure 2 shows the phasor diagrams of the salient-pole synchronous generator and the salient-pole synchronous motor.In Figure 2, φ, λ, and ψ0 represent the power factor angle, the electric angle between the stator phase current and the d-axis, and the internal power factor angle, respectively.U is the stator voltage, Ė is the air-gap electromotive force, İ is the stator current, and Xσ is the armature leakage reactance.
The stator resistance value of the large salient-pole synchronous alternator is small (usually in the order of 0.001 Ω) and can be ignored.Under the premise of ignoring the stator resistance, the equation for the phase voltage is as follows: Before calculating the electromagnetic field of the machine, the first step is to determine the field current If and the power factor angle φ.The iterative process of determining If and φ is shown in Figure 3. Firstly, pre-select a reasonable If and λ, which are substituted into the electromagnetic model for the first calculation step.Then, the waveform of center air-gap flux density B can be calculated from the finite element analysis.The amplitude of the air-gap flux density Bm and the angle between air-gap magnetic flux Φ and d-axis can be obtained using the fast Fourier transform of B. As Φ leads Ė by 90°, the angle ψ1 between Ė and the q-axis equal to the angle between Φ and the d-axis.

Electromagnetic Field Calculation
The pumped storage generator motor has two working conditions: generating and pumping, which can be treated as a salient-pole synchronous generator and a salient-pole synchronous motor, respectively.Figure 2 shows the phasor diagrams of the salient-pole synchronous generator and the salient-pole synchronous motor.In Figure 2, ϕ, λ, and ψ 0 represent the power factor angle, the electric angle between the stator phase current and the d-axis, and the internal power factor angle, respectively.
. U is the stator voltage, Ė is the air-gap electromotive force, İ is the stator current, and X σ is the armature leakage reactance.

2023, 11, x FOR PEER REVIEW 5 of 21
According to the phasor relationship in Figure 2, the angles β and α can be derived as follows: where E = 4.44fNkWΦ, Φ = 2Bmτlef/π, lef is the effective axial length of the stator, N is the number of turns in series for each phase of the stator, τ is the pole pitch, and kW is the fundamental winding factor.The angle φ can be calculated according to ψ1, α, and the pre-selected λ.The stator voltage U can be derived from the phasor diagram as follows: In this study, Xσ is calculated by the analytical method.The calculated value of Xσ is 0.1183 for the generating operation and 0.1340 for the pumping operation.Afterward, de- The stator resistance value of the large salient-pole synchronous alternator is small (usually in the order of 0.001 Ω) and can be ignored.Under the premise of ignoring the stator resistance, the equation for the phase voltage is as follows: Before calculating the electromagnetic field of the machine, the first step is to determine the field current I f and the power factor angle ϕ.The iterative process of determining I f and ϕ is shown in Figure 3 where UN is the stator phase voltage under rated load, φN is the rated power factor, and ε is the convergence tolerance, usually taken as 0.5% or less.
After the iterative process, the electromagnetic field distributions of the machine under rated generating operation and pumping operation are shown in Figure 4.The magnetic field distribution of the machine under generating operation is similar to that of pumping operation.The maximum magnetic flux density occurs at the edge position of the rotor pole shoe, with a value of 2.18 T. The magnetic flux density of the pumping operation is slightly higher than that of the generating operation.According to the phasor relationship in Figure 2, the angles β and α can be derived as follows: where E = 4.44fNk W Φ, Φ = 2B m τl ef /π, l ef is the effective axial length of the stator, N is the number of turns in series for each phase of the stator, τ is the pole pitch, and k W is the fundamental winding factor.The angle ϕ can be calculated according to ψ 1 , α, and the pre-selected λ.The stator voltage U can be derived from the phasor diagram as follows: In this study, X σ is calculated by the analytical method.The calculated value of X σ is 0.1183 for the generating operation and 0.1340 for the pumping operation.Afterward, determine whether the computing residuals of U and ϕ meet the precision tolerance of Equation (5).If U and ϕ meet the precision requirements, the iterative computation stops, obtaining I f and ϕ.If U and ϕ do not meet accuracy requirements, adjust the I f and λ by I f ' = k 1 I f and λ' = λ + k 2 (ϕ − ϕ N ) and replace the original value with the corrected value in the electromagnetic model to re-calculate U and ϕ.Repeat the above steps until U and ϕ meet the accuracy requirements of Equation (5).
where U N is the stator phase voltage under rated load, ϕ N is the rated power factor, and ε is the convergence tolerance, usually taken as 0.5% or less.
After the iterative process, the electromagnetic field distributions of the machine under rated generating operation and pumping operation are shown in Figure 4.The magnetic field distribution of the machine under generating operation is similar to that of pumping operation.The maximum magnetic flux density occurs at the edge position of the rotor pole shoe, with a value of 2.18 T. The magnetic flux density of the pumping operation is slightly higher than that of the generating operation.

Rotor Losses Calculation
The electromagnetic losses of the rotor region consist of three parts: copper loss of the field winding, additional losses on the rotor pole surface, and eddy current loss of the damping winding.The copper loss of the field winding is computed by I 2 R and the calculated result is 286.2 kW.

Eddy Current Loss of Damping Bar
Due to the open slot of the stator, the uneven distribution of permeance along the circumferential direction of the stator has an impact on the air-gap main magnetic field, making a series of harmonic magnetic fields produced in the air-gap.Then, the relative motion is generated between those fields and the rotor, inducing an eddy current within the damping winding circuit.The magnetic fields produced by such induced current can offset the effect of additional magnetic fields, but the eddy current loss is also generated in damping winding as a result.
The eddy current distribution of the machine's rated generating operation and rated pumping operation is obtained by the 2D transient electromagnetic field, as shown in Figure 5.The maximum current density is centralized at the top of the damping bar.The eddy current loop also exists within a single damping bar.The eddy current loss induced in damping bars can be calculated according to the eddy current density: where n e is the element amount of one damping bar, J e is the eddy current density of one single element (in A/m 2 ), ∆ e is the element area (in m 2 ), L pf is the axial length of the damping bar (in m), σ is conductivity (in S/m).
(a) (b) Figure 6 shows the calculated results of the eddy current loss of the damping bar.In under rated generating operation, the minimum eddy current loss value of the 1st damping bar closed to the windward side is 132 W, while the maximum eddy current loss value of the 8th damping bar closed to the leeward side is 492 W.Under rated pumping operation, the minimum eddy current loss value of the 1st damping bar closed to the leeward side is 171 W, while the maximum eddy current loss value of the 8th damping bar closed to the windward side is 582 W. In addition, the eddy current loss from the 2nd to the 7th damping bar is close under two working conditions.The armature reaction has an impact on the electromagnetic field distribution, making each damping bar generate different eddy current loss.The total eddy current loss of the damping bar under rated generating operation and the rated pumping operation is 2775 W and 2877 W, respectively.The overall loss of the damping winding is relatively small under two working conditions.

Additional Loss of Rotor Pole Surface
The additional loss induced in the rotor pole surface includes three parts [22]:

Additional Loss of Rotor Pole Surface
The additional loss induced in the rotor pole surface includes three parts [22]: 1.
Additional loss of rotor pole surface under no-load rated voltage P FeP (in Watt): where K δ1 is the air-gap coefficient of the stator tooth, B δ is the average flux density of the air-gap (in Tesla), t 1 is the pitch of the stator tooth (in cm), Z is the stator slot number, n N is the synchronous speed (in r/min), p is the number of pole pair, A p is the surface area of the rotor pole (in cm), and K 1 and K 2 are correction factors which are related to material of the rotor pole or the damper winding if the rotor pole is composed of solid forged steel, the value of K 1 is 23.3; if the rotor pole is composed of solid cast steel, the value of K 1 is 17.5; if the rotor pole is composed of sheets of silicon steel, the value of K 1 is determined in accordance with Figure 7.The value of K 2 can be chosen in accordance with Table 2.

Additional Loss of Rotor Pole Surface
The additional loss induced in the rotor pole surface includes three parts [22]: 1. Additional loss of rotor pole surface under no-load rated voltage PFeP (in Watt): where Kδ1 is the air-gap coefficient of the stator tooth, Bδ is the average flux density of the air-gap (in Tesla), t1 is the pitch of the stator tooth (in cm), Z is the stator slot number, nN is the synchronous speed (in r/min), p is the number of pole pair, Ap′ is the surface area of the rotor pole (in cm), and K1 and K2 are correction factors which are related to material of the rotor pole or the damper winding if the rotor pole is composed of solid forged steel, the value of K1 is 23.3; if the rotor pole is composed of solid cast steel, the value of K1 is 17.5; if the rotor pole is composed of sheets of silicon steel, the value of K1 is determined in accordance with Figure 7.The value of K2 can be chosen in accordance with Table 2.

2.
Additional loss caused by harmonic magnetomotive force (MMF) of the stator windings P kν (in Watt): where p is the number of pole pair, n B , l B , and A B are the number, length (in cm), and crosssectional area (in cm 2 ) of the damping winding, respectively, C B is the ratio of resistivity of damping winding and cooper, δ is the length of the air-gap (in cm), K δ is Carter's coefficient, λ δ = (t 2 − b sh )/K δ δ , b sh is the opening width of the damping winding (in cm), t 2 is the damping winding pitch (in cm), λ S = 0.5(1 − 0.64τ/l E β), β is the ratio of the short pitch, and d B is the damping winding diameter (in cm), F δ = 1.6δK δ B δ .Additional loss caused by harmonic MMF of the stator teeth P 2νk (in Watt): where p is the number of pole-pair, k is a correction factor, K δ is Carter's coefficient.
Machines 2023, 11, 901 The calculated additional loss of the rotor pole surface is summarized in Table 3.The P Fep of the machine is the same under both generating operation and pumping operation and accounts for the largest proportion of total additional losses.P kν and P 2νk account for a relatively small proportion of the total additional loss, while proportions of P kν and P 2νk are slightly larger under pumping operation.The proportions of P kν and P 2νk under the generating operation are 3.2% and 20.2%, respectively.The proportions of P kν and P 2νk under the pumping operation are 3.8% and 24.3%, respectively.In order to investigate the velocity distribution of cooling air and the temperature distribution of the rotor, the 3D fluid-thermal coupled model is built in this section.Considering the ventilation system structure of the real machine, as shown in Figure 8a, the 3D fluid-thermal coupled model of the machine rotor is established as shown in Figure 8b.Cooling air from the heat exchanger flows to both ends of the stator end-windings to firstly cool the stator end-windings (path ①-②-③).And then through the rotor end region, the cooling air flows into the rotor supporter under the effect of rotor rotation.Afterward, the cooling air flows into the rotor yoke radial ventilation ducts and then enters the air region between two adjacent rotor poles to cool the field winding (path ③-④-⑤ ).At last, cooling air flows into the stator radial ventilation ducts through the air-gap, taking away the heat generation of the armature winding and stator core, and then re-enters the heat exchanger to dissipate heat (path ⑤-⑥-①).
According to the symmetry of the heat transfer path, which includes the upper ventilation channel and lower ventilation channel, the model of one rotor pole pair and half length of the machine is established  Cooling air from the heat exchanger flows to both ends of the stator end-windings to firstly cool the stator end-windings (path 1 -2 -3 ).And then through the rotor end region, the cooling air flows into the rotor supporter under the effect of rotor rotation.Afterward, the cooling air flows into the rotor yoke radial ventilation ducts and then enters the air region between two adjacent rotor poles to cool the field winding (path 3 -4 -5 ).At last, cooling air flows into the stator radial ventilation ducts through the air-gap, taking away the heat generation of the armature winding and stator core, and then re-enters the heat exchanger to dissipate heat (path 5 -6 -1 ).
According to the symmetry of the heat transfer path, which includes the upper ventilation channel and lower ventilation channel, the model of one rotor pole pair and half length of the machine is established.The surface of the rotor supporter is treated as the inlet of the solved rotor model, and cooling air inlets of the stator ventilation ducts are treated as the outlets of the established rotor model.The heat transfer path 4 -5 -6 is modeled in the established model.Twenty-eight rotor outlets, i.e., the inlets of the stator ventilation ducts, and one outlet of the rotor end region are established in the rotor fluid-thermal coupled model.Based on the fluid dynamic and heat transfer theories, the solved region satisfies the mass conservation equation, momentum equation, and energy conservation equation.The governing equations of the fluid-thermal coupled model of the pumped storage generator motor are given as follows [23][24][25].
where ρ is fluid density (in kg/m 3 ), t is the time (in s), u is the velocity vector of the fluid, u i , u j , and u k are the x y z axial components of u (in m/s), λ cond is the thermal conductivity (in W/(m•K)), P is the fluid pressure (in Pa), µ is the turbulent viscosity coefficient (in kg/m•s), T is the temperature (in • C), c is the specific heat (in J/(kg•K)), S T is the internal heat source of the fluid (in W/m 3 ), and S u , S v , and S w are the heat sources of the fluid.
The following assumptions are proposed for the 3D fluid-thermal coupled field analysis: 1.
The Reynolds number is widely used as a criterion to determine whether the fluid flow state is obviously greater than 2300.The fluid in the ventilation system of the machine is in a turbulent state, and the standard k-ε model is used to calculate the flow of cooling air.The Boltzmann equations with kinetic energy of turbulence k and diffusion factor ε are as follows: where ε is the diffusion factor (in m 2 /s 3 ).G k is the generation rate of the turbulence, µ t is the turbulent viscosity coefficient, C 1ε and C 2ε are constants, and σ k and σ ε are Planck constants of turbulent flow.

2.
The influences of aerostatics buoyancy and gravity on fluid motion are neglected; 3.
The flow velocity of cooling air is far less than the acoustic velocity, and thus cooling air is regarded as an incompressible fluid.
The boundary conditions of the 3D fluid-thermal coupled model are as follows: 1. S1, the inlet of the rotor supporter, is defined as the velocity inlet boundary condition.
The outlet velocities of the four heat exchangers of the machine are measured on the real machine, as listed in  The surfaces (in red as shown in Figure 8b) at both sides of the rotor model are defined as periodic boundary conditions; 4.
Considering the symmetry in an axial direction, the axially central section (in blue as shown in Figure 8b) is considered to be an adiabatic surface; 5.
The room temperature is considered as 13 • C.
The additional loss of the rotor pole surface and eddy current loss of damping winding obtained from the electromagnetic calculation are transferred to 3D fluid-thermal coupled field analysis as heat generation sources.The cooling air velocity distribution and the temperature distribution of rotor components are determined after 135 iterations using the commercial software package Fluent.

Analysis of Cooling Air Velocity Distribution
In the fully air-cooled large alternator, the heat dissipation of the crucial heat generation components, such as field winding, primarily depends on the cooling effect of cooling air between two adjacent rotor poles.To gain a better understanding of the cooling air velocity distribution between two adjacent rotor poles, the velocity distribution along the axial direction is studied, as illustrated in Figure 9.After cooling air enters the rotor through the rotor supporter, a portion of cooling air enters the rotor pole space through the radial ventilation ducts of the rotor yoke, and another portion enters the rotor pole space through the end zone.Total outlets in the rotor-solved region include twenty-eight surfaces of stator radial ventilation ducts and one end region outlet, as shown in Figure 9.The cooling air from rotor ventilation ducts flows approximately radial to the 1st~28th outlets.However, the flow of the cooling air from the end region is also relatively large, and it will be forced to continue moving along the axial direction until it reaches the cooling air outlet and gradually flows into the stator.
The cooling air velocity at the ventilation duct of the rotor yoke is 65 m/s (as marked in white label in Figure 9), while the cooling air velocity gradient at the end area is relatively large, which ranges from 35 m/s to 50 m/s.As a result of the Coriolis force, the cooling air velocity between adjacent rotor poles gradually grows along the radial direction, and the cooling air velocity of the end region grows from 65 m/s to approximately 90 m/s.  Figure 10 gives the cooling air velocity vector distribution of different radial sections between two adjacent rotor poles at 28 cooling air outlets of the rotor.Along the axial direction, the 1st radial section (corresponding to the 1st outlet) is located near the axial center of the rotor, while the 28th radial section (corresponding to the 28th outlet) is at the rotor end region.The cooling air velocity between two adjacent rotor poles exhibits an increasing trend along the axial direction from the rotor center to the end region.However, it should be noted that the maximum cooling air velocity is 99.35 m/s which is located in the 26th radial section.The studied machine rotor has a total of nine yoke ventilation ducts along the axial direction.In this paper, the adiabatic boundary is applied at the axial center, and, hence, 4.5 rotor yoke ventilation ducts are established in the rotor model.Along the axial direction, the first ventilation duct is located near the end region, while the fifth ventilation duct is located at the axial center of the rotor.Figure 11 shows the average and maximum cooling air velocity of the ventilation channel of the rotor.The average air speeds of five ventilation ducts range from 63.3 to 70.5 m/s, while the maximum airspeeds range from 99.7 m/s to 103.0 m/s.The differences between the average air speeds and maximum airspeeds are 1.2 m/s and 3.3 m/s, respectively.Therefore, the velocity distribution of cooling air in each ventilation channel is very even to improve the heat dissipation condition of the central part of the field winding.In the first radial section, it is observed that the velocity direction of cooling air between two adjacent rotor poles closes to the circumferential direction instead of the radial direction.This is owing to the joint influence of rotor rotation and large rotor diameter.Compared to the leeward side, the cooling airflow direction at the windward side is closer to the Y direction (radial direction).At those locations with a high velocity of cooling air, the flow direction of cooling air shifts to the radial direction.The velocity distribution pattern of cooling air at the 26th radial section is similar to that of the 1st radial section.However, due to the higher cooling air velocity of the 26th radial section, the effect of the Coriolis force is more significant, resulting in higher cooling air velocity at an upper position between two adjacent rotor poles.
The studied machine rotor has a total of nine yoke ventilation ducts along the axial direction.In this paper, the adiabatic boundary is applied at the axial center, and, hence, 4.5 rotor yoke ventilation ducts are established in the rotor model.Along the axial direction, the first ventilation duct is located near the end region, while the fifth ventilation duct is located at the axial center of the rotor.Figure 11 shows the average and maximum cooling air velocity of the ventilation channel of the rotor.The average air speeds of five ventilation ducts range from 63.3 to 70.5 m/s, while the maximum airspeeds range from 99.7 m/s to 103.0 m/s.The differences between the average air speeds and maximum airspeeds are 1.2 m/s and 3.3 m/s, respectively.Therefore, the velocity distribution of cooling air in each ventilation channel is very even to improve the heat dissipation condition of the central part of the field winding.

Analysis and Experiment Validation of the Rotor Thermal Field
Figure 12 shows the temperature distribution of various rotor components under rated generating operation and rated pumping operation.As can be seen in Figure 12, the temperature distribution law of the rotor is as follows: along the axial direction, the rotor temperature gradually decreases from the axial center to the end region; along the circumferential direction, the temperature gradually decreases from the leeward side to the windward side; along the direction, the temperature gradually decreases from the rotor shaft to the rotor pole surface.
As per Figure 8, it can be found that after entering the rotor supporter, a large portion of the cooling air enters the rotor pole space through the end region, which will apparently help to cool the end part of the field winding.This part of cooling air will continue to move along the axial direction from the end region to the axial center region, and gradually flow out the rotor to take away the heat generation.The heat dissipation condition in the end region is better than that of the axial center region.Therefore, the rotor temperature gradually decreases from the axial center to the end region.
Because of the rotor rotation, the airflow velocity and corresponding heat transfer coefficient at the windward side of the rotor pole are apparently higher than that of the leeward side.Therefore, all rotor components, including field winding and damper bars, etc., exhibit temperature variation law that the temperature gradually decreases from the leeward side to the windward side along the circumferential direction.
Compared to the rotor part which is cooled by the flowing air in the rotor radial ventilation ducts, the heat dissipation condition of the rotor pole surface is much better, which is cooled by the cooling air in the air-gap and has a larger effective heat transfer area.Therefore, along the radial direction, the temperature gradually decreases from the rotor shaft to the rotor pole surface.Figure 13 shows the maximum temperature, minimum temperature, and average temperature of various rotor parts.It can be found from Figure 13a that in rated generating operation, the high-temperature area of the rotor is predominantly centralized on the field

Analysis and Experiment Validation of the Rotor Thermal Field
Figure 12 shows the temperature distribution of various rotor components under rated generating operation and rated pumping operation.As can be seen in Figure 12, the temperature distribution law of the rotor is as follows: along the axial direction, the rotor temperature gradually decreases from the axial center to the end region; along the circumferential direction, the temperature gradually decreases from the leeward side to the windward side; along the radial direction, the temperature gradually decreases from the rotor shaft to the rotor pole surface.

Analysis and Experiment Validation of the Rotor Thermal Field
Figure 12 shows the temperature distribution of various rotor components under rated generating operation and rated pumping operation.As can be seen in Figure 12, the temperature distribution law of the rotor is as follows: along the axial direction, the rotor temperature gradually decreases from the axial center to the end region; along the circumferential direction, the temperature gradually decreases from the leeward side to the windward side; along the radial direction, the temperature gradually decreases from the rotor shaft to the rotor pole surface.
As per Figure 8, it can be found that after entering the rotor supporter, a large portion of the cooling air enters the rotor pole space through the end region, which will apparently help to cool the end part of the field winding.This part of cooling air will continue to move along the axial direction from the end region to the axial center region, and gradually flow out the rotor to take away the heat generation.The heat dissipation condition in the end region is better than that of the axial center region.Therefore, the rotor temperature gradually decreases from the axial center to the end region.
Because of the rotor rotation, the airflow velocity and corresponding heat transfer coefficient at the windward side of the rotor pole are apparently higher than that of the leeward side.Therefore, all rotor components, including field winding and damper bars, etc., exhibit temperature variation law that the temperature gradually decreases from the leeward side to the windward side along the circumferential direction.
Compared to the rotor part which is cooled by the flowing air in the rotor radial ventilation ducts, the heat dissipation condition of the rotor pole surface is much better, which is cooled by the cooling air in the air-gap and has a larger effective heat transfer area.Therefore, along the radial direction, the temperature gradually decreases from the rotor shaft to the rotor pole surface.Figure 13 shows the maximum temperature, minimum temperature, and average temperature of various rotor parts.It can be found from Figure 13a that in rated generating operation, the high-temperature area of the rotor is predominantly centralized on the field As per Figure 8, it can be found that after entering the rotor supporter, a large portion of the cooling air enters the rotor pole space through the end region, which will apparently help to cool the end part of the field winding.This part of cooling air will continue to move along the axial direction from the end region to the axial center region, and gradually flow out the rotor to take away the heat generation.The heat dissipation condition in the end region is better than that of the axial center region.Therefore, the rotor temperature gradually decreases from the axial center to the end region.
Because of the rotor rotation, the airflow velocity and corresponding heat transfer coefficient at the windward side of the rotor pole are apparently higher than that of the leeward side.Therefore, all rotor components, including field winding and damper bars, etc., exhibit temperature variation law that the temperature gradually decreases from the leeward side to the windward side along the circumferential direction.
Compared to the rotor part which is cooled by the flowing air in the rotor radial ventilation ducts, the heat dissipation condition of the rotor pole surface is much better, which is cooled by the cooling air in the air-gap and has a larger effective heat transfer area.Therefore, along the radial direction, the temperature gradually decreases from the rotor shaft to the rotor pole surface.
Figure 13 shows the maximum temperature, minimum temperature, and average temperature of various rotor parts.It can be found from Figure 13a that in rated generating operation, the high-temperature area of the rotor is predominantly centralized on the field winding.The maximum temperature, minimum temperature, and average temperature of the field winding are 75.6 • C, 39.3 • C, and 57.6 • C, respectively.
winding.The maximum temperature, minimum temperature, and average temperature of the field winding are 75.6 °C, 39.3 °C, and 57.6 °C, respectively.From a comparison of Figure 12a,b, it can be observed that the temperature distribution law of various rotor components under rated pumping operation is the same as that under generating operation.Therefore, the temperature distribution under pumping operation will not be further analyzed in the following analysis of temperature distribution.From a comparison of Figure 13a,b, it is notable that the overall average temperature rise of the rotor under rated pumping operation is 3 °C higher than that under generating operation.The reason is that the copper loss of the field winding and the additional loss are greater under pumping operation.
Figure 14 illustrates the temperature variation of the field winding.Owing to the effect of the high cooling air velocity in the end region of the field winding, the temperature of the field winding end region is relatively low, and the temperature range is 39.3~56 °C.To investigate the temperature variation rule of the field winding center, one radial sample plane is selected at the field winding center on both sides, respectively.The black radial sample plane 1 is located inside the rotor field winding at the windward side, while the red radial sample plane 2 is located inside the rotor field winding at the leeward side.In addition, 12 sample lines numbered 1~12 are built from the lower part to the upper part of the field winding along the radial direction, as illustrated in Figure 14a.From a comparison of Figure 12a,b, it can be observed that the temperature distribution law of various rotor components under rated pumping operation is the same as that under generating operation.Therefore, the temperature distribution under pumping operation will not be further analyzed in the following analysis of temperature distribution.From a comparison of Figure 13a,b, it is notable that the overall average temperature rise of the rotor under rated pumping operation is 3 • C higher than that under generating operation.The reason is that the copper loss of the field winding and the additional loss are greater under pumping operation.
Figure 14 illustrates the temperature variation of the field winding.Owing to the effect of the high cooling air velocity in the end region of the field winding, the temperature of the field winding end region is relatively low, and the temperature range is 39.3~56 • C. To investigate the temperature variation rule of the field winding center, one radial sample plane is selected at the field winding center on both sides, respectively.The black radial sample plane 1 is located inside the rotor field winding at the windward side, while the red radial sample plane 2 is located inside the rotor field winding at the leeward side.In addition, 12 sample lines numbered 1~12 are built from the lower part to the upper part of the field winding along the radial direction, as illustrated in Figure 14a.
The temperature variation of 12 sample lines in both the windward side and the leeward side of the field winding are shown in Figure 14b,c.It can be seen that the temperature of the field winding center gradually decreases from the bottom to the top.As a result of the high cooling airflow velocity of the windward side, the temperature of the windward side is significantly lower than the leeward side at the identical position.Furthermore, the temperature variation range of the field winding on the leeward side is larger, which ranges from 50 • C to 75 • C. The temperature distribution on the windward side is relatively uniform, which ranges from 50 • C to 60 • C.
The temperature distribution of the rotor pole surface Is shown In Figure 15.The additional loss of one single rotor pole shoe surface under rated generating operation is 3.3 kW, and the rotor pole surface is in good heat dissipation condition.Therefore, the surface temperature of each rotor pole is relatively low, with a maximum temperature of 43 • C and a minimum temperature of 22 • C. The temperature distribution of the rotor pole surface gradually decreases from the leeward side to the windward side.Additionally, the temperature difference of the windward and leeward sides gradually decreases as it approaches the rotor end region.The temperature difference near the rotor center is 8 • C, while near the end region, the temperature difference reduces to 2 • C. The temperature variation of 12 sample lines in both the windward side and the leeward side of the field winding are shown in Figure 14b,c.It can be seen that the temperature of the field winding center gradually decreases from the bottom to the top.As a result of the high cooling airflow velocity of the windward side, the temperature of the windward side is significantly lower than the leeward side at the identical position.Furthermore, the temperature variation range of the field winding on the leeward side is larger, which ranges from 50 °C to 75 °C.The temperature distribution on the windward side is relatively uniform, which ranges from 50 °C to 60 °C.
The temperature distribution of the rotor pole surface is shown in Figure 15.The additional loss of one single rotor pole shoe surface under rated generating operation is 3.3 kW, and the rotor pole surface is in good heat dissipation condition.Therefore, the surface temperature of each rotor pole is relatively low, with a maximum temperature of 43 °C and a minimum temperature of 22 °C.The temperature distribution of the rotor pole surface gradually decreases from the leeward side to the windward side.Additionally, the temperature difference of the windward and leeward sides gradually decreases as it approaches the rotor end region.The temperature difference near the rotor center is 8 °C, while near the end region, the temperature difference reduces to 2 °C. Figure 16a,b shows the temperature distribution of eight damping bars in one rotor pole and the temperature variation curve in the central axis position of each damping bar under the generating operation, respectively.In Figure 16b, the horizontal axis refers to the axial distance, while zero refers to the axial center position.The temperature of the damping bar decreases from the axial center to the end in a large range of 24.2~40.6°C.The maximum axial temperature difference of one single damping bar is up to 15.0 °C.In the circumferential direction, the additional loss of the 8th damping bar is the maximum.As the eighth damping bar is located in the leeward side of the machine's rotation direction and the cooling effect is poorer, the maximum temperature of the eighth damping bar is the highest at 40.8 °C.The additional loss of the first damping bar is the minimum.As   Figure 16a,b shows the temperature distribution of eight damping bars in one rotor pole and the temperature variation curve in the central axis position of each damping bar under the generating operation, respectively.In Figure 16b, the horizontal axis refers to the axial distance, while zero refers to the axial center position.The temperature of the damping bar decreases from the axial center to the end in a large range of 24.2~40.6°C.The maximum axial temperature difference of one single damping bar is up to 15.0 °C.In the circumferential direction, the additional loss of the 8th damping bar is the maximum.As the eighth damping bar is located in the leeward side of the machine's rotation direction and the cooling effect is poorer, the maximum temperature of the eighth damping bar is the highest at 40.8 °C.The additional loss of the first damping bar is the minimum.As the first damping bar is located the closest to the windward side, its maximum temperature is the lowest among the eight damping bars, which is 37.0 °C.The temperature of the eight damping bars decreases from the eighth to the first damping bar.The temperature rise test is conducted at the pump storage power plant of the studied machine in this paper to validate the accuracy of the calculated results.The resistance test approach is widely employed to measure the average temperature rise of the field winding for large alternators.The average temperature rise can be obtained based on the relationship that the winding resistance varies correspondingly with the temperature rise.The average temperature rise could be determined as: The temperature rise test is conducted at the pump storage power plant of the studied machine in this paper to validate the accuracy of the calculated results.The resistance test approach is widely employed to measure the average temperature rise of the field winding for large alternators.The average temperature rise can be obtained based on the relationship that the winding resistance varies correspondingly with the temperature rise.The average temperature rise could be determined as: where R 1 is the winding resistance under the ambient temperature (in Ω), R 2 is the field winding resistance under the measured temperature (in Ω), t 1 is the ambient temperature corresponding to R 1 , t 2 is the cooling medium temperature corresponding to R 2 , and K is a constant, for copper winding it is 235.The operating parameters of the machine under the temperature test are summarized in Table 5.The electromagnetic loss under the testing condition is re-calculated through the electromagnetic model and the 3D fluid-thermal coupled model is also re-solved with associated losses.The calculated average temperature of the field winding under the testing condition is 46.1 • C, and the error is 1.1% compared to the measured value of 45.6 • C, thus verifying the calculation accuracy in this article.

Conclusions
This paper calculates the additional losses of the rotor pole surface and the eddy current loss of the different damping bars based on the established 2D transient electromagnetic model.Furthermore, the 3D coupled fluid-thermal model is established to analyze the complex airflow velocity distribution of cooling air and the temperature distribution of various rotor components.The consistency between the calculated temperature of the field winding and the test result verifies the accuracy of the model.The research results of this paper are expected to give empirical insights into the heat dissipation effect of such a full air-cooled ventilation system and help to develop it to be used in the future design of a larger capacity hydrogenerator.The following conclusions can be obtained.

1.
The additional loss of the rotor pole surface and the eddy current loss of the damping winding account for a relatively small proportion of the total rotor loss.The additional loss and the eddy current loss are slightly greater in rated pumping operation compared to the rated generating operation.2.
In the space between two adjacent rotor poles, the cooling air velocity gradually decreases along the axial direction from the end region to the rotor center.Simultaneously, the cooling air velocity gradually increases along the radial direction from the rotor shaft to the rotor pole surface, and the maximum cooling air velocity is 99.35 m/s.The cooling air velocity distribution of five ventilation ducts of the rotor is uniform, and the air speeds difference of the average and maximum are 1.2 m/s and 3.3 m/s, respectively.

3.
The temperature distribution law of various rotor components of the pumped storage generator motor is as follows: the temperature of the rotor area gradually decreases from the axial center to the end region along the axial direction; the temperature gradually decreases from the leeward side to the windward side along the circumferential direction; and the temperature gradually decreases along the radial direction from the rotor shaft to the rotor pole surface.The highest temperature of the rotor appears in the field winding zone.Specifically, in rated generating operation and rated pumping operation, the highest temperatures are 75.6 • C and 78.3 • C, respectively.Uniform cooling air distribution can effectively improve the heat dissipation performance of the rotor area.
. Firstly, pre-select a reasonable I f and λ, which are substituted into the electromagnetic model for the first calculation step.Then, the waveform of center air-gap flux density B can be calculated from the finite element analysis.The amplitude of the air-gap flux density B m and the angle between air-gap magnetic flux .Φ and d-axis can be obtained using the fast Fourier transform of B. As .Φ leads Ė by 90 • , the angle ψ 1 between Ė and the q-axis equal to the angle between .Φ and the d-axis.

Figure 3 .
Figure 3. Determination of I f and ϕ.Symbol ∆ is the algebraic deviation between iterative and rated values, and k 1 and k 2 are convergence factors.0 < k 1 < 1 and 0 < k 2 < 1.

Figure 4 .
Figure 4. Magnetic flux density distributions of the machine.(a) Rated generating operation; (b) rated pumping operation.Figure 4. Magnetic flux density distributions of the machine.(a) Rated generating operation; (b) rated pumping operation.

Figure 4 .
Figure 4. Magnetic flux density distributions of the machine.(a) Rated generating operation; (b) rated pumping operation.Figure 4. Magnetic flux density distributions of the machine.(a) Rated generating operation; (b) rated pumping operation.

Figure 6 .
Figure 6.Eddy current loss of damping bar.(a) Rated generating operation; (b) rated pumping operation.

Figure 6 .
Figure 6.Eddy current loss of damping bar.(a) Rated generating operation; (b) rated pumping operation.
. The surface of the rotor supporter is treated as the inlet of the solved rotor model, and cooling air inlets of the stator ventilation ducts are treated as the outlets of the established rotor model.The heat transfer path ④-⑤-⑥ is modeled in the established model.Twenty-eight rotor outlets, i.e., the inlets of the stator ventilation ducts, and one outlet of the rotor end region are established in the rotor fluidthermal coupled model.Based on the fluid dynamic and heat transfer theories, the solved

Machines 2023 , 21 Figure 9 .
Figure 9. Cooling air velocity distribution on axial section between two adjacent rotor poles.Figure 9. Cooling air velocity distribution on axial section between two adjacent rotor poles.

Figure 9 .
Figure 9. Cooling air velocity distribution on axial section between two adjacent rotor poles.Figure 9. Cooling air velocity distribution on axial section between two adjacent rotor poles.

Figure 9 .
Figure 9. Cooling air velocity distribution on axial section between two adjacent rotor poles.

Figure 10 .
Figure 10.Cooling air velocity vector distribution on radial sections between two adjacent rotor poles.

Figure 10 .
Figure 10.Cooling air velocity vector distribution on radial sections between two adjacent rotor poles.

Figure 11 .
Figure 11.Cooling air velocity at the outlets of the rotor yoke ventilation ducts.

Figure 12 .
Figure 12.Temperature distribution of various rotor components.(a) Rated generating operation; (b) rated pumping operation.

Figure 11 .
Figure 11.Cooling air velocity at the outlets of the rotor yoke ventilation ducts.

Figure 11 .
Figure 11.Cooling air velocity at the outlets of the rotor yoke ventilation ducts.

Figure 12 .
Figure 12.Temperature distribution of various rotor components.(a) Rated generating operation; (b) rated pumping operation.

Figure 12 .
Figure 12.Temperature distribution of various rotor components.(a) Rated generating operation; (b) rated pumping operation.

Figure 13 .
Figure 13.Maximum, minimum, and average temperature of various rotor components.(a) Rated generating operation; (b) rated pumping operation.

Figure 13 .
Figure 13.Maximum, minimum, and average temperature of various rotor components.(a) Rated generating operation; (b) rated pumping operation.

Figure 14 .
Figure 14.Temperature variation of the field winding.(a) Temperature distribution; (b) temperature variation at windward side; (c) temperature variation at leeward side.

Figure 14 . 21 Figure 15 .
Figure 14.Temperature variation of the field winding.(a) Temperature distribution; (b) temperature variation at windward side; (c) temperature variation at leeward side.Machines 2023, 11, x FOR PEER REVIEW 18 of 21

Figure 15 .
Figure 15.Temperature distribution of the rotor pole surface under rated generating operation.

Figure
Figure 16a,b shows the temperature distribution of eight damping bars in one rotor pole and the temperature variation curve in the central axis position of each damping bar under the generating operation, respectively.In Figure 16b, the horizontal axis refers to the axial distance, while zero refers to the axial center position.The temperature of the damping bar decreases from the axial center to the end in a large range of 24.2~40.6• C. The maximum axial temperature difference of one single damping bar is up to 15.0 • C. In the

Figure 15 .
Figure 15.Temperature distribution of the rotor pole surface under rated generating operation.

Table 1 .
Geometric Dimensions and Electromagnetic Parameters.

Table 3 .
Additional loss of rotor pole surface.

Calculation and Analysis of Rotor Fluid-Thermal Coupled Field
4.1.Establishment of Rotor Fluid-Thermal Model

Table 4 .
The average outlet velocity is 2.83 m/s and the outlet area of each heat exchanger is 3.198 m 2 .Thus, the total circulation flow of cooling air is calculated as Q = 2.83 m/s × 3.198 m 2 × 8 = 72.4m 3 /s (the total amount of the heat exchanger is eight).The inlet velocity of S1 is calculated by Q/S, and the value is 4.68 m/s, where S is the surface area of the model inlet;

Table 5 .
Operating parameters of the machine under the test.