Coupling calculation and analysis of three-dimensional temperature and fluid field for High-power High-Speed Permanent Magnet Machine

: In order to accurately estimate the temperature rise for high-power high-speed permanent magnet machines (HSPMMs), a novel temperature calculation method based on multi-physics co-simulation is proposed. Based on the computational fluid dynamics (CFD) and heat transfer theory, the computation model of fluid-solid-heat coupling heat transfer is established, and the coupled field is calculated using finite volume method with fundamental assumptions and corresponding boundary conditions. With the influences from temperature gradient and water flow rate considered, the heat transfer coefficients of water pipe surfaces is obtained by the application of the inverse iteration method. Thus, HSPMM temperature and fluid field can be simulated numerically by the finite volume methods, while the spatial temperature distributions for the machine main components are analyzed in this paper. The 1.12 MW, 18000 r/min HSPMM is prototyped with experiments conducted on it, while the test data are then compared with the calculated results, which validate the correctness of the solution method of the coupled field.


Introduction
High speed permanent magnet machines (HSPMMs) are developing towards the high power density characteristic. HSPMMs demonstrate their advantages in weight and volume aspects when compared with the conventional machines at the same power rate. HSPMMs are also capable to directly drive the high speed load, which greatly benefits the reduction of weight, volume and noise for the whole drive system. However, temperature rise is normally considered as a major concern for the application of HSPMMs [1]- [3]. Firstly, both temporal and spatial harmonics act in HSPMM air gap, which cause large machine iron losses, copper losses and rotor eddy losses. It is also noted that the air frictional losses by the rotor high speed operation are also found significantly increased when compared with such losses of the conventional machines [4]. Moreover, HSPMMs are more vulnerable to suffer from rotor over-heating problems due to the compact size of HSPMMs. In addition, carbon fiber sleeve is necessary to be wrapped outside rotor PMs for mechanical considerations [5]，which further limits rotor heat dissipation capability. Therefore, high power HSPMMs are characterized with strong-coupling, multi-constraints and nonlinear properties, while their temperature field calculations are more complicated than the conventional machines.
The main methods for temperature rise calculation can be concluded as lumped parameter thermal path method and temperature field method [6]- [7]. Many scholars have analyzed and calculated fluid fields and temperature fields for various machines. A 3-D temperature field calculation method for PM motors based on finite formula method and numerical heat transfer theory is developed in [8], thus the calculation time can be effectively saved. The methods for determining the heat transfer coefficient of the winding end is investigated in [9], while the influences of air duct structure deformation on fluid field and temperature field are also discussed. Literature [10] analyzes the fluid field in the air-cooled system of a megawatt HSPMM. However, in the simulation model of temperature field, the windings and insulation layers in stator slots are difficult to be modeled as the actual structure utilized in the machine [11]. The surface heat dissipation coefficient of liquid-cooled for high-density PM machine is decided based on the fluid flow state, the fluid physical properties and the geometric parameters of the heat dissipation surface. In addition, these parameters also restrict each other. The heat dissipation coefficient is not a constant, and it is difficult to determine accurately. At present, empirical methods are usually used to approximate the key parameters such as thermal conductivity and convection heat transfer coefficient [12], and there is still a gap towards the accurate calculation of machine temperature rise.
Based on the characteristics of high-power HSPMM, a novel temperature rise calculation process based on multi-physics co-simulation is firstly proposed in this paper. The temperature rise distributions for machine components are also obtained. A 1.12 MW, 18000 rpm HSPMM is prototyped and tested at rated condition, while the temperature test results agree well with the simulation results. Then, the fluid-solid-heat coupling and inverse iteration method are developed to calibrate the heat transfer coefficients for the prototype in order to perform advanced thermal analyses and achieve the desirable cooling performance for the machine.

Temperature calculation method based on multiphysics co-simulation
The design parameters of the 1.12MW, 18000r/min HSPMM are shown in Table1. The rotor topology adopts the surface-mounted PM structure, the PM in each pole is divided into three segments and pole clearances are filled with pole filler. Rotor sleeve is made of carbon fiber composite material. The combined stator water-cooled with rotor air-cooled are also designed for the machine.
The calculation method of HSPMM's loss is different from an ordinary motor. The loss analysis and calculation results of this prototype have been discussed in detail in [13] and [14]. The conductive properties of materials vary with temperatures, which may affect their performances in electrical machines. Therefore, in this calculation and analysis process, the influences of working temperature on the electromagnetic and thermal properties of materials in the HSPMM are fully considered, which is performed through the temperature iteration.
Calculation steps of temperature rise for high-power HSPMM based on multi-physical field analysis can be concluded as: 1) assume the working temperature of the machine and the convection heat transfer coefficient on the surface of the waterway, and calculate the machine losses; 2) calculate the temperature distribution for the machine based on the fluid-solid-heat coupling method, and determine the convection heat transfer coefficients by the inverse iterative method; 3) if the calculated temperature is inconsistent with the assumed temperature, the iterative calculation is then repeated until the requirements are met. The temperature judgment could be written as: Where Tk is assuming the machine working temperature; Tk ' is the calculated. The temperature iteration calculation flow is shown in Fig. 1.

Solving domain physical models and boundary conditions
Due to the circumferential symmetrical structure of the machine, a range of stator slot distance of the HSPMM is selected for the 3-D solution domain model, while the winding ends are straightened in the model. The inner conductors are equivalent to uniform copper strips，while the insulation layer and air are equivalent to insulating entities. The solving domain model is presented in Fig. 2. The turbulence model is used to solve the flow field in the machine for the stable operation of HSPMM. The fluid in the motor is considered to be incompressible. The boundary conditions for the solution of the temperature field and the fluid field can be regarded as:  Inlet wind speed is 14m/s. Outlet pressure is 1 standard atmospheric pressure.
 The convection heat transfer coefficient need to be given for the water channel surface.
Inlet water temperature is assumed to be 20°C, and the flow rate is 3m 3 /h. The rotor rotational speed is 180 m/s.

Mathematical model
According to computational fluid dynamics (CFD) and numerical heat transfer, fluid flow follows the laws of physical conservation [15]: Mass conservation equation is given as: where ρ is fluid density, t is time, V is velocity vector.
Momentum conservation equation is expressed as: where F is surface force, P is fluid pressure, η is viscosity coefficient of the fluid.
Energy conservation equation is written as: where λ is the thermal conductivity and is heat flux.

Determination of coefficient of thermal conductivity and convective heat transfer
The HSPMM's losses are time-varying, and the coefficient of thermal conductivity of machine material is also affected by the variation of the environment and temperature. Therefore, machine thermal performance estimation requires the combination of numerical calculation and finite element analysis method. According to the principle of equivalent heat source, the silicon steel sheet of stator core is completely pressed, and the core material is a thermally conductive anisotropic material. The thermal coefficients of radial and tangential thermal paths are 42W/(mC), and the axial thermal conductivity is 0.45W/(mC) [16]; The carbon fiber sleeve is Waterway anisotropic in thermal properties, with a radial and tangential thermal coefficient of 0.7 W/(mC) and an axial thermal conductivity of 4.8 W/(mC) [17].
The heat transfer between the internal cooling water and the waterway surface is performed as forcing convection. According to the fluid similarity theory, the convective heat transfer coefficient of the waterway surface is expressed as [17]: To sum up: Where  is the convective heat transfer coefficient of the waterway surface; Nu is Nusselt number which indicates the strength of convective heat transfer capacity; d is the equivalent diameter; λ is the fluid thermal conductivity coefficient; Re is the Reynolds number which indicates fluid flow characteristics; Pr is the Prandtl number at normal temperature to reflect the relative magnitude of fluid momentum diffusion capacity and heat diffusion capacity; k is the coefficient considering flow character and geometric parameters of heat dissipation surface; ρ is the density of the fluid, v is the fluid flow rate, and μ is the fluid viscosity.
Accurate calculation of the convective heat transfer coefficient is critical to the correctness of the thermal analysis results. The convective heat transfer coefficient of the waterway surface is usually a constant value in previous calculation. In reality, the convective heat transfer coefficient along the axial direction of the waterway surface is different with the changes of motor temperature and water flow temperature. This paper proposes an iterative correction method to determine the convective heat transfer coefficient of the waterway surface. The assumed convective heat transfer coefficient along the axial direction is given as 1   6  . (The 6 data points along the waterway are shown in Fig.3).
As the water flow moves along the axial direction, the temperature gradually increases and the convective heat transfer coefficient decreases. According to (7) that Fourier heat conduction law and Newton's law of heat dissipation, the appropriate convective heat transfer coefficient is solved by multiple iteration correction.
where T is the node temperature of the casing surface, Tf is the fluid temperature.
As shown in equation (7), nonlinear relationship is found between the convective heat transfer coefficient and the temperature. Therefore, the convective heat transfer coefficient of the channel along the axial direction is approximated as: where T1 and T2 are the temperature of different nodes on waterway surface, ΔX is the distance between two nodes, ' i  is the approximate value of the convective heat transfer coefficient of the i-th node, αi is the correction value of the convective heat transfer coefficient of the i-th node, β is the relaxation factor that is iterative for the convective heat transfer coefficient.
According to the assumed convective heat transfer coefficient, the temperature of each node can be obtained by fluid-solid-heat coupling method. According to (8)  ), and all the calculations will be performed once again until the convergence conditions in (9) are satisfied. Namely, when the temperature difference between two iterations of the same node is less than 0.5%, and the difference of the convective heat transfer coefficient is also less than 0.5%, the iteration ends. The judgment of iteration calculation for convective heat transfer coefficient can be expressed as: where TS is the temperature measured on the real machine in each point highlighted in Fig.3, and TW is the corresponding calculated temperature.
After the above iterative correction, the curve of the convective heat transfer coefficient along the axial direction is fitted under the experimental conditions when the water flow velocity (v0) is 0.5m/s. The curves of the convective heat transfer coefficient along axial direction at different water flow velocities are obtained by (6) and presented in Fig. 4.

Analysis of temperature field calculation results
After getting the temperature field and the fluid field coupled result of the physical model under 25°C ambient temperature, the temperature rise distribution in the whole machine can be obtained as presented in Fig. 5. It is learned from Fig.5 that the highest temperature rise of the stator is 80°C near the end of the winding outlet, and the rotor maximum temperature rise is 145°C and locates in rotor middle region. In this case, PMs may suffer from irreversible demagnetization for machine at long-term operation. The radial temperature rise distributions are shown in Fig.6. Rotor carbon fiber sleeve has poor thermal conductivity and anisotropic thermal property as its radial thermal conductivity is much smaller than the axial thermal conductivity. Therefore, local temperature rise of the PM is found quite high, due to the difficulty of transmitting the large amount of heat generated by the rotor loss in the radial direction. Moreover, there is a risk of PM demagnetization, which affects the running performance and reliability of the machine.
The rotor temperature rise distribution along the axial Temperature rise (C) direction is shown in Fig.7. The rotor temperature rising near the air inlet and outlet are 87°C and 107°C, respectively. Fig.8 presents the relationship between the convective heat transfer coefficient of the waterway and the temperature rise of the case. As the temperature increases, the convective heat transfer coefficient of the waterway decreases, and the heat dissipation capacity decreases.

Temperature rise test
The temperature rise experiment for the prototype is carried out. The test platform is shown in Fig.9. The machine is powered by variable frequency drive (VFD) with voltage. It drives an asynchronous generator with the coupling between the two machines achieved using a gear box. Hence, the high speed PMSM output can be reflected in the power generated by the asynchronous generator. The VFD is a nine level high voltage inverter with vector control method applied to achieve control in both machine speed and power during high speed operation. The nine level high voltage inverter has the advantage of powering the machine with low harmonics, and it is a desirable solution for high power rate high speed machine drive. Fig.10 presents the heat flux density distribution of each component and the position where temperature sensor locates. The rotor temperature is measured by an infrared temperature sensor. The prototype machine is tested under the rated conditions (current and speed). Fig.11 presents the experimental data and the winding phase current of the machine at the rated speed (18000 rpm) and rated load condition. The calculated ambient temperature is the same as actual one in the temperature rise test. When the machine temperature rise is less than 0.1C per hour, temperature reaches a steady state. The measured values compared to the proposed method and the lumped parameter thermal networks method in [6]. The compared results are shown in Table.3. Through comparative analysis, it is observed that the calculation temperature rise results of the rotor and stator are basically consistent with the measured values. Therefore, the correctness of the numerical simulation results is verified. The temperature field calculation method proposed in this paper is more

Synchronous generator Prototype
Air cooling Water cooling Fig.9 Prototype temperature rise test platform

Effect of axial ventilation speed on rotor cooling
Effective rotor cooling is greatly desired due to the high heat flux density of the high-power HSPMM (see Fig.10) and the temperature rising. According to the analysis of the fluid field, the axial ventilation speed determines the rotor heat dissipation. However, excessive wind speed also causes extra air friction loss and noise. In order to investigate the influences of ventilation speeds on rotor cooling, different ventilation speeds have been considered from 1m/s to 50m/s. The calculation results are given in Fig.12.
When the wind speed varies in the range from 1 m/s to 10 m/s, the temperature rise of the rotor drops sharply. However, after the wind speed exceeding 10m/s, the cooling effect is little improved by the wind speed increasing, and the heat dissipation capacity is close to saturation. Therefore, both machine thermal performance and maintenance cost can benefit from the selection of the optimum ventilation speed for HSPMM. Fig. 13 shows the relationship between the velocity of rotor surface and the ventilation speed. The velocity VR of the rotor surface is approximately linear with the optimum ventilation velocity Vv, which can be expressed as following:

Influence of internal air duct height and air gap variation on optimal wind speed
In this section, a sensitivity analysis is performed to investigate the impact of air gap and air duct on the optimal ventilation speed. When other conditions keep unchanged, only the size of air gap and the height of the inner duct are changed, the optimal ventilation speed Vv, which corresponding to rotor velocity VR, is calculated by the fluid field and temperature field coupling method. It is also compared with the calculated result of analytical method by (9). Fig.14 and Fig.15 present the optimal wind speed variation with respect of air duct clearance and duct height, respectively, when both fields coupling method and analytical method are utilized.
It is found that the change of the air gap and the air duct height have almost no influence on the optimal matching relationship between the rotor velocity VR and the optimal ventilation speed Vv, while the rotor surface velocity is the decisive factor. The simulated values by fields coupling method match well with the calculated values by analytical method, which proving (9) is widely applicable.

Conclusion
In this paper the temperature calculation method based on multi-physics co-simulation is proposed. The problem of temperature accurate calculation is effectively solved due to interaction between loss and temperature rise, and the following conclusions are drawn: The coupled calculation model of the 3-D temperature field and the fluid field for the high-power HSPMM is established. The temperature field and fluid field analysis of the machine are carried out, and the correctness of the calculation model is verified by the temperature rise experimental test of the prototype, which providing a reference for the high-power HSPMM cooling system design.
Considering the nonlinear relationship among cooling water flow rate, temperature gradient and convective heat transfer coefficient, the reverse convection iteration method is proposed and applied to determine the convective heat transfer coefficient of the waterway surface, and the accuracy of the temperature field calculation is improved.
Increasing the ventilation speed will improve the cooling effect for the rotor under forced ventilation condition. However, the heat dissipation capacity is saturated when the wind speed reaching a certain value. There is an optimum ventilation speed for a rotor surface velocity, which is not affected by the air gap and the air duct height. A reasonable ventilation speed should be selected to improve the cooling effect for HSPMM. Inner air duct height 36mm Inner air duct height 39mm