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

Abstract

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 significant 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 field and fluid–thermal coupled field of a pumped storage generator motor are analyzed. The 2D transient time-step finite element model of the electromagnetic field of a pumped storage generator motor is established, and the eddy current loss of damping bars of the rotor is calculated by the finite element method. The additional loss of the rotor pole surface is calculated by analytical method. The mathematical and geometric models of the 3D fluid–thermal coupled field of the pumped storage generator motor are established and calculated. The complex fluid 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 field 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.

1. 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 high-quality 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 effect is not as good as the direct cooling method, such as the water-cooled or hydrogen-cooled method, etc. Cooling air inside the machine is completely driven by the rotor’s own rotation to take away the heat generation of the machine. Cooling air exhibits a complicated distribution within the alternator which is cooled by the fully air-cooled method. The temperature rise of various components of the machine is directly influenced by the distribution and flow rule of cooling air. Therefore, it is of great significance to study the fluid field distribution of cooling air and the thermal field of the machine rotor.

In recent years, some researchers have studied the cooling mechanism and temperature distribution of the large salient-pole synchronous alternator [4,5,6,7,8,9,10,11,12,13]. Bersch et al. developed a novel thermal and electromagnetic coupled model to obtain optimal ventilation design. To reduce the peak temperature of the stator winding, the radial ventilation of the stator was optimized and designed [14]. Han and Li et al. conducted a study on the temperature distribution and fluid distribution of the end region of a 330-MW water–hydrogen-cooled turbogenerator using the transient electromagnetic field and the fluid–solid conjugated heat transfer method. In order to reduce the temperature rise of the machine end region, various methods such as changing the air gap spacer height, the shelter board width, and the electric shield material were studied [15,16,17]. Wang et al. established a 3D fluid–thermal coupled model of synchronous condenser end region and analyzed the effect of the cooling medium velocity on the temperature of various components of the rotor [18]. Fan et al. established a 3D thermal field model of the hydrogenerator and studied the influence of the structural parameters on the damper winding heat generation of the hydrogenerator. The study found that the damper bar pitch can affect the losses and heat generation of damper bars [19,20]. Carounagarane et al. analyzed the thermal behavior of a large hydrogenerator under continuous overload operation. The result shows that the temperature rise of the generator with Class F insulation is acceptable under 20% continuous overload operation [21]. However, the thermal analysis is carried out by applying the empirical radiative and convective heat transfer coefficients as the boundary conditions in the thermal model. The temperature calculation result of this method heavily depends on the accuracy of the heat transfer coefficients obtained from empirical formulas, which is often challenging in large hydrogenerators with complex structures.

Numerous scholars have conducted research on the cooling and heat dissipation problems of large synchronous machines. However, the research on complex fluid field distribution and thermal field of a large synchronous machine that employs a full air-cooled fanless cooling system is not yet comprehensive. The cooling air in such a ventilation system is driven by the rotation effect of the rotor support and rotor pole. The cooling fan is canceled to significantly reduce the ventilation loss and the machine enjoys the operation economy and stability. Cooling air is forced to move within such a large dimension rotating rotor, which apparently exhibits a complex fluid distribution within the rotor. To gain a better understanding of the heat dissipation effect of such a ventilation system, and develop it to be used in larger capacity hydrogenerator, a comprehensive investigation on heat generation, cooling airflow distribution, and temperature field of the machine with a fully air-cooled fanless cooling system is highly desired. The coupled fluid–thermal analysis is the most accurate and effective approach to study this issue.

This paper takes a fully air-cooled 145 MW (generating)–170 MW (pumping) pumped storage generator motor as an example, and the 2D transient electromagnetic field and 3D fluid–thermal coupled field of the pumped storage generator motor are calculated and analyzed. Firstly, the 2D time-step finite element model of transient electromagnetic field based on the real machine is established, and the eddy current loss of the damping bar is obtained by the finite element method. The additional loss on the rotor pole surface is calculated using an analytical method. Then, the 3D fluid–thermal coupled model of the fully air-cooled pumped storage generator motor is established, and corresponding governing equations and boundary conditions of the solving region are given. The calculated losses are applied as heat generation sources in the 3D fluid–thermal coupled model. Finally, the flow rule of cooling air and the temperature distribution of the machine rotor are analyzed in detail.

2. Model Establishment and Calculation of Transient Electromagnetic Filed

2.1. 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. Geometric Dimensions and Electromagnetic Parameters.

Parameter	Generating Operation	Pumping Operation
Power rating (MW)	145	170
Rated voltage (kV)	13.8	13.8
Rated current (kA)	6.9	7.82
Power factor	0.88	0.91
Efficiency (%)	98.4	98.5
Frequency (Hz)	50
Rated speed (r/min)	200
Air gap (mm)	28
Number of stator slot	360
Number of pole pair	15
Stator outer diameter (mm)	8100
Stator inner diameter (mm)	7240
Width of stator slot (mm)	24.5
Height of stator slot (mm)	159.5
Axial length (mm)	2050
Effective stator axial length (mm)	1611
Rotor outer diameter (mm)	7184
Rotor pole inner diameter (mm)	6644
Width of field winding (mm)	90
Thickness of field winding (mm)	8
Turns of field winding	27
Stator core weight (kg)	109,700
Stator copper weight (kg)	18,186
Rotor copper weight (kg)	27,104

.

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:

• 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;
• 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:

$$\{\begin{array}{l}\frac{\partial }{\partial x}(\frac{1}{\mu }\cdot \frac{\partial {\mathit{A}}_z}{\partial x})+\frac{\partial }{\partial y}(\frac{1}{\mu }\cdot \frac{\partial {\mathit{A}}_z}{\partial y})=-{\mathit{J}}_z+\sigma \frac{\partial {\mathit{A}}_z}{\partial t}\hfill \\ {\mathit{A}}_z|AB={\mathit{A}}_z|CD=0\hfill \\ {\mathit{A}}_z|AD={\mathit{A}}_z|BC\hfill \end{array}$$

(1)

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.

2.2. 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 ψ₀ represent the power factor angle, the electric angle between the stator phase current and the d-axis, and the internal power factor angle, respectively. $\dot{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:

$$\dot{U}=\dot{E}-j\dot{I}{X}_{\sigma }$$

(2)

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. 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 $\dot{\Phi }$ and d-axis can be obtained using the fast Fourier transform of B. As $\dot{\Phi }$ leads Ė by 90°, the angle ψ₁ between Ė and the q-axis equal to the angle between $\dot{\Phi }$ and the d-axis.

According to the phasor relationship in Figure 2, the angles β and α can be derived as follows:

$$\{\begin{array}{l}\beta ={\psi }_1+\lambda \\ \alpha =\arctan \frac{I{X}_{\sigma }\sin \beta }{E-I{X}_{\sigma }\sin \beta }\end{array}$$

(3)

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 ψ₁, α, and the pre-selected λ. The stator voltage U can be derived from the phasor diagram as follows:

$$U=\sqrt{E^2+{(I{\mathrm{X}}_{\sigma })}^2-2EI{\mathrm{X}}_{\sigma }\mathrm{c}\mathrm{o}\mathrm{s}\beta }$$

(4)

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₁I_f and λ’ = λ + k₂ (φ − φ_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).

$$\{\begin{array}{l}\left|\frac{U-U_N}{U_N}\right|\times 100\%<\epsilon \\ \left|\frac{\cos \phi -\cos {\phi }_N}{\cos {\phi }_N}\right|\times 100\%<\epsilon \end{array}$$

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

3. 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²R and the calculated result is 286.2 kW.

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

$$P_i={\displaystyle \sum _{e=1}^{n_e}{P}_e=}{\displaystyle \sum _{e=1}^{n_e}{J}_e^2{\Delta }_e}{L}_{pf}{\sigma }^{-1}$$

(6)

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²), Δ_e is the element area (in m²), L_pf is the axial length of the damping bar (in m), σ is conductivity (in S/m).

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.

3.2. 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 under no-load rated voltage P_FeP (in Watt):

$$P_{Fep}=\left(K_1+K_2\right){\left[\frac{\left(K_{\mathsf{\delta }1}-1\right)B_{\mathsf{\delta }}{t}_1}{{10}^3}\right]}^2{\left[\frac{Z{n}_N}{{10}^4}\right]}^{1.5}2p{A}_p{}^{\prime }\times {10}^{-4}$$

(7)

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₁ 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₁ and K₂ 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₁ is 23.3; if the rotor pole is composed of solid cast steel, the value of K₁ is 17.5; if the rotor pole is composed of sheets of silicon steel, the value of K₁ is determined in accordance with Figure 7. The value of K₂ can be chosen in accordance with

Table 2. Determination of K₂.

Damper Winding	K2
Without damper winding	0
Red copper	2.5
Brass	10

.

2.

Additional loss caused by harmonic magnetomotive force (MMF) of the stator windings P_kν (in Watt):

$$P_{Kv}=\frac{3.7p{n}_B{l}_B{C}_B{K}_{f300}}{A_B}{\left(\frac{\tau }{\delta K_{\delta }}\right)}^2\times {\left(\frac{K_{\delta }{x}_{ad}(F_{\delta }/2000)}{({\lambda }_{\delta }/1-\cos {\delta }_B)+2{\lambda }_S}\right)}^2\times {10}^{-2}$$

(8)

where p is the number of pole pair, n_B, l_B, and A_B are the number, length (in cm), and cross-sectional area (in cm²) 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, ${\lambda }_{\delta }=(t_2-b_{sh})/K_{\delta }{\delta }^{\prime }$, b_sh is the opening width of the damping winding (in cm), t₂ is the damping winding pitch (in cm), ${\lambda }_S=0.5\left(1-0.64\tau /l_E\beta \right)$, β is the ratio of the short pitch, $K_{f300}=1.96d_B\sqrt{1/C_B}$, and d_B is the damping winding diameter (in cm), $F_{\delta }=1.6\delta K_{\delta }{B}_{\delta }$.

Additional loss caused by harmonic MMF of the stator teeth P_2νk (in Watt):

$$P_{2\nu k}=k^{\prime }{\left[\frac{2p{X}_{ad}}{Z_1(K_{\delta }-1)}\right]}^2{P}_{Fep}\times {10}^3$$

(9)

where p is the number of pole-pair, k′ is a correction factor, K_δ is Carter’s coefficient.

The calculated additional loss of the rotor pole surface is summarized in

Table 3. Additional loss of rotor pole surface.

Condition	PFeP/kW	Pkν/kW	P2νk/kW	Total Additional Loss/kW
Generating	75.73	3.14	19.96	98.83
Pumping	75.73	4.03	25.60	105.36

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

4. Calculation and Analysis of Rotor Fluid–Thermal Coupled Field

4.1. Establishment of Rotor Fluid–Thermal Model

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

$$\{\begin{array}{l}\frac{\partial \rho }{\partial t}+div\left(\rho \mathit{u}\right)=0\\ \frac{\partial \left(\rho u_i\right)}{\partial t}+div\left(\rho \mathit{u}{u}_i\right)=div\left(\mu \cdot grad{u}_i\right)-\frac{\partial P}{\partial x}+S_u\\ \frac{\partial \left(\rho u_j\right)}{\partial t}+div\left(\rho \mathit{u}{u}_j\right)=div\left(\mu \cdot grad{u}_j\right)-\frac{\partial P}{\partial y}+S_v\\ \frac{\partial \left(\rho u_k\right)}{\partial t}+div\left(\rho \mathit{u}{u}_k\right)=div\left(\mu \cdot grad{u}_k\right)-\frac{\partial P}{\partial z}+S_w\\ \frac{\partial \left(\rho T\right)}{\partial t}+div\left(\rho \mathit{u}T\right)=div\left(\frac{{\lambda }_{\mathrm{cond}}}{c}\cdot gradT\right)+S_T\end{array}$$

(10)

where ρ is fluid density (in kg/m³), 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³), 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:

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

$$\{\begin{array}{l}\frac{\partial }{\partial t}(\rho k)+\mathrm{div}(\rho k\mathit{u})=\mathrm{div}[(\mu +\frac{{\mu }_t}{{\sigma }_k})\mathrm{grad}k]+G_k-\rho \epsilon \\ \frac{\partial }{\partial t}(\rho \epsilon )+\mathrm{div}(\rho \epsilon \mathit{u})=\mathrm{div}[(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }})\mathrm{grad}\epsilon ]+C_{1\epsilon }\frac{\epsilon }{k}{G}_k-\rho C_{2\epsilon }\frac{{\epsilon }^2}{k}\end{array}$$

(11)

where ε is the diffusion factor (in m²/s³). 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:

• 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

  Table 4. Ventilation measurement results.

  Number of Heat Exchangers	Outlet Velocity (m/s)	Average Velocity (m/s)	Outlet Area (m2)
  1	2.78	2.83	3.198
  2	3.01
  7	2.76
  8	2.73

  . The average outlet velocity is 2.83 m/s and the outlet area of each heat exchanger is 3.198 m². Thus, the total circulation flow of cooling air is calculated as Q = 2.83 m/s × 3.198 m² × 8 = 72.4 m³/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;

2.

S2, the outlet surface of the rotor, is defined as a natural outflow boundary condition;

3.

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.

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

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.

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

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 8^th 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:

$$\theta =\frac{R_2-R_1}{R_1}(K+t_1)+t_1-t_2$$

(12)

where R₁ is the winding resistance under the ambient temperature (in Ω), R₂ is the field winding resistance under the measured temperature (in Ω), t₁ is the ambient temperature corresponding to R₁, t₂ is the cooling medium temperature corresponding to R₂, 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. Operating parameters of the machine under the test.

Parameter	Value
Stator voltage (kV)	13.96
Stator current (kA)	6.11
Power factor	0.993
Field winding voltage (kV)	0.149
Field winding current (kA)	1.316

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

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

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