Abstract
With the development of electrical machines for achieving higher performance and smaller size, heat generation in electrical machines has also increased. Consequently, the temperature rise in electrical machines causes unexpected heating of components and makes it difficult to operate properly. Therefore, in the development of electrical machines, the accurate evaluation of temperature increase is important. In the thermal design of electrical machines, heat-conduction analysis using the heat-transfer boundary set on the surface of a heated target has been frequently performed. However, because the heat-transfer coefficient is dependent on various factors, it is often determined based on experimental or numerical simulation results. Therefore, setting the heat-transfer coefficient to a constant value for the surface of the heated target degrades the analysis accuracy because the actual phenomenon cannot be modeled. To enhance the accuracy of the heat-transfer coefficient, the coupled electromagnetic field with heat-conduction analysis finite element method (FEM), thermal-fluid analysis using FEM, and the highly simplified marker and cell method is applied to the estimation of the distribution of the heat-transfer coefficient. Moreover, to accurately calculate the localized heat-transfer coefficient, the temperature distribution and flow velocity distribution around the heated target are analyzed in the induction-heating apparatus.
1 Introduction
In the development of electrical machines, components are placed with high density for high performance and miniaturization. Consequently, the heat generation density rises in electrical machines, causing unexpected heating of components, which makes it difficult to properly operate. Thus, in the development of high-quality and smaller-sized electrical machines, evaluating the temperature-rise characteristics in the actual driving state of electrical machines is important.
As a simulation technique to evaluate the temperature rise in the actual driving state of electrical machines, coupled electromagnetic-field and heat-conduction analysis using the heat-transfer boundary set on the surface of the heated target has been frequently carried out [1,2]. However, because the heat-transfer coefficient changes depending on location owing to various factors (e.g., shape, material of the heated target, and temperature of the surroundings), it is often determined such that the discrepancy between the measurement data and numerical result is minimized [3]. Therefore, from the viewpoint of advanced prediction of machine performance, the advantage of numerical simulation is not fully maximized. Furthermore, it is a rare case that the heat-transfer coefficient is uniformly distributed on the surface of the heated target, and modeling the local change in the heat-transfer coefficient from actual measurement results is difficult.
Therefore, thermal fluid analysis is often adopted as a technique to model the local heat-transfer effect. Several numerical simulation techniques have been proposed to simulate fluid phenomena. In the particle method [4], which is often adopted in liquid analysis, the Lagrange coordinate system [5,6] is applied; thus, the physical quantity of only the position where the particle exists is calculated. The Euler coordinate system is often applied in cases in which the physical quantity of the entire analysis area is required, e.g., analysis of the flow of gas. In the Euler coordinate system, the marker and cell (MAC) [7], [8] method, which separately solves the flow velocity and pressure fields, is frequently used. In the fluid analysis, the highly simplified MAC (HSMAC) [9] method, which reduces the matrix calculation (e.g., matrix-vector product) by applying the Newton method to the pressure calculation stage in the MAC method, is adopted. HSMAC method is more advantageous than other fluid analysis methods. When the HSMAC method is compared with the finite element method (FEM)-based thermal fluid analysis, the difficulty of implementation of FEM is more difficult than that of the HSMAC method. Furthermore, the analysis stability in FEM is strongly dependent on the combination of discretization order for the pressure and velocity [10]. Therefore, the HSMAC method is applied to the solver for the thermal-fluid field in this paper.
In this study, the behavior of natural convection derived from the induction heating phenomenon is simulated by the coupled electromagnetic-field with heat conduction and thermal-fluid analysis using FEM and HSMAC method. Furthermore, the evaluation method of the local heat-transfer coefficient [11] using the HSMAC method is demonstrated by comparison with actual measurements and numerical results.
2 Analysis method
2.1 Eddy-current analysis using edge-based FEM
The weak form of the edge-based FEM using magnetic vector potential A and electric scalar potential ϕ is defined as follows:
where Ωall, Ωc, and Ωe are the whole analysis region, input current region, and eddy current region, respectively. The material constants ν and σ are the magnetic reluctivity and electric conductivity, respectively. The physical quantity J 0 is the magnetizing current density, N i and N j are the edge-based and node-based shape functions, respectively. The index i and j are the computational points defined for the edges and nodes, respectively. In Ωc, J 0 is defined as −grad ψ where ψ is the electric scalar potential defined in the magnetizing winding region and ϕ is calculated to satisfy div(−grad ϕ) = 0. The eddy current is formulated using −(∂ A /∂t + grad ϕ). The tangential component of A is continuous between the next elements and ϕ is defined only in Ωc. The convergence characteristic of the ICCG method in A − ϕ method is better than that of A method [12].
2.2 Heat-conduction analysis using node-based FEM
The weak form of the heat-conduction equation is defined as follows:
where i is the index of computation node, Ωh is the thermal conduction area, λ is the thermal conductivity, ρ is the density, c is the specific heat capacity, T is the temperature, Q is the heat-source density, and Γ S is the surface of the heated target. In the solid region to be analyzed for heat-conduction, heat exchange occurs in the fluid because of the heat-transfer phenomenon in (3). Thus, the boundary integral term in (3) is changed to the following:
where h is the heat-transfer coefficient, and T f is the fluid temperature.
2.3 Thermal-fluid analysis using HSMAC method
The Navier–Stokes equation, which is the governing equation for incompressible fluid, is defined as (5). The external force term in the third term on the right is applied to the Boussinesq approximation. Next, the continuity equation is defined as (6). The energy equation is described as follows(7):
where u is the velocity, p is the pressure, η is the viscosity coefficient, β is the coefficient of volumetric expansion, and g is the gravitational acceleration.
In this study, the HSMAC method is adopted as the thermal-fluid analysis technique. The finite difference method (FDM) is applied to the spatial discretization in the HSMAC method, and an explicit method is applied to the discretization in time marching. To discretize the advection term, the first-order upwind difference [13] is applied. In addition, a staggered grid is adopted as a computational mesh.
A staggered grid evaluates the scalar and vector quantities at different computational points. A computational point is defined at the center of a cell for scalar quantities such as pressure p and temperature T. Furthermore, the computational points of vector quantities, such as the velocity in the x-direction u x , velocity in the y-direction u y , and velocity in the z-direction u z are defined as the position shifted by a half-cell from the center of a cell. As an example, a 2D staggered grid is shown in Figure 1. The discretized distance in each direction is shown in Figure 2, where i and j are the index of computational points in the x- and y-directions, respectively, and Δx f, Δy f, and Δz f and Δx b, Δy b, and Δz b are the distances from the computational points to the discrete points that go forward and backward in the x-, y-, and z-directions, respectively.
Staggered grid.
Discretized distance. (a) Whole view. (b) x–y plane. (c) z–x plane.
For the discretization in (7), because thermal conductivity λ cannot be treated as a constant in a multimedia interface, the diffusion term is discretized as follows:
where k is the index of computational points in the z-direction. The continuity in the normal direction of the heat flux at the boundary between the heated target and air is imposed by (8).
2.4 Calculation of heat transfer coefficient
Fourier’s law is defined as (9):
where q is the heat flux. Newton’s cooling law is defined as (10):
T w is the surface temperature and T f is the fluid temperature. The fluid temperature T f refers to the vicinity of the point at which the heat-transfer coefficient is calculated. Because a direction of the heat flux in the boundary is an outward normal direction, vector q can be a scalar q. By eliminating q from (9) and (10), the local heat-transfer coefficient is formulated as follows:
2.5 Weakly coupled analysis
The flowchart of the coupled electromagnetic-field and thermal-fluid analysis is shown in Figure 3. Initial conditions are set (step 1). The time step of the electromagnetic-field analysis is advanced (step 2). The eddy-current loss obtained from the electromagnetic-field analysis is substituted for Q in (7) (steps 3 and 4). The time step of the thermal-fluid analysis is advanced (step 5). The temporary velocity and pressure correction values are calculated, and the temporary velocity is corrected, which are the computational steps of the HSMAC method (steps 6–8). The HSMAC method is iterated until the continuity equation is satisfied (step 9). The temperature distribution is calculated using the known velocity obtained from the HSMAC method (step 10). The conductivity is set to the constant value in the fluid calculation loop (steps 5–10). Since the thermal-fluid analysis is calculated in the low-temperature range, conductivity of the heat target can be defined as linear characteristic. The discretized equations using the explicit method in the time direction are parallelized using OpenMP in steps 6 and 10. The temporary velocity is calculated in step 6, and the temporary temperature is calculated in step 10. Thermal-fluid analysis is carried out until the time of the thermal-fluid analysis is larger than that of the electromagnetic-field analysis (step 11). The analysis of the proposed method is repeated until the analysis time of the thermal fluid reaches the maximum time (step 12).
Flowchart of the weakly coupled electromagnetic-field and thermal-fluid analysis.
3 Analysis model and experimental instrument
3.1 Induction-heating model
As shown in Figure 4, the SS400 is heated by eddy current induced from the induction-heating coil. The electromagnetic-field analysis is carried out in one-half of the model of Figure 4. The temperature is numerically evaluated at point α and γ near the surface of the SS400. Expanded polyethylene (EPE) is used as a heat-insulating material inside the coil. The adiabatic boundary condition is applied outside the air region.
Induction-heating model. (a) x–y plane, (b) whole view.
In the electromagnetic-field analysis, discretization with skin depth taken into account is required. Therefore, to prevent the increase in the computational scale, the eddy current loss density is interpolated and applied to the heat-conduction analysis. As an example, one element on the conductor surface is shown in Figure 5. The interpolation formula is defined as (12). The eddy current loss density calculated by electromagnetic field analysis is interpolated by (12), and the heat-source density of each element is obtained. The eddy current loss of each element is calculated by (13). The electromagnetic field is formulated in the time domain in equation (13). To analyze the eddy current distribution in Ωe, the width of one finite element in the direction of skin depth is set to less than 1/3 of skin depth. However, when the mesh to its fine mesh is applied to the thermal-fluid analysis, the calculation time will be unrealistic. Therefore, the heat source density is interpolated to the coarse mesh which is applied to the thermal-fluid analysis.
where Q e is the averaged heat-source density of the element obtained by interpolation, n is the number of elements included in the range to be interpolated, P i is the eddy-current loss of each element, and Δx, Δy, and Δz are the lengths of the element obtained by interpolation in the x-, y-, and z-directions, respectively.
Interpolation of eddy-current loss.
3.2 Measuring instruments for induction heating
Figure 6 shows the experimental apparatus. IH-1M, manufactured by NAVIO Corp., is used as the induction-heating apparatus. CT6863-05, manufactured by HIOKI Corp., is used as the current sensor. TBS1042, manufactured by Tektronix Corp., is used as the oscilloscope. FLIR ONE for iOS Gen 3 PRO, manufactured by FLIR Corp., is used for the thermography. The current waveform used in the electromagnetic-field analysis is a sine wave obtained by averaging the waveform observed by the current sensor in the oscilloscope. The amplitude of the sine wave is 28.8 A, and the frequency is 28.5 kHz.
Experimental instrument. (a) Whole view, (b) induction heating coil.
4 Analysis results
4.1 Analysis condition
Table 1 lists the types of coupled analysis techniques. The abbreviation “FE–FE” is defined as the coupled method in which both the electromagnetic field and heat conduction are solved by FEM. In FE–FE, heat-conduction analysis without considering the heat-transfer effect is applied. “FE–FD” is defined as the coupled method in which the electromagnetic field is solved by FEM, and heat conduction is solved by FDM. In FE–FD, heat-conduction analysis without considering the heat-transfer effect is applied. “FE–FD–HS” is defined as the coupled method in which the electromagnetic field is solved by FEM, heat conduction is solved by FDM, and thermal-fluid analysis is solved by HSMAC method. All software is our own codes. Table 2 lists the computational scale of each analysis, convergence criterion of the preconditioned CG method, and Newton method. Table 3 lists the material constants. For the numerical simulation, a computer is used with the following specifications: CPU, Intel Core i7-6850 (six cores, 3.6 GHz), and 128 GB RAM.
Combinations for coupled analysis
| Coupled method | EM | Thermal | Fluid |
|---|---|---|---|
| FE–FE | FEM | FEM | |
| FE–FD | FEM | FDM | |
| FE–FD–HS | FEM | FDM | HSMAC |
Computational scale and convergence criterion
| EM | Heat conduction | HSMAC method | |||||
|---|---|---|---|---|---|---|---|
| u x | u y | u z | p | T | |||
| DoF | 7,92,162 | 2,28,528 | 1,71,166 | 1,69,432 | 1,68,656 | 1,75,616 | 2,28,528 |
| Nonzero | 1,44,95,918 | Without coefficient matrix | |||||
| ε CG | 10−6 | ||||||
| ε NR | 10−3 | 10−3 | |||||
| Δt | 1.096 × 10−6 | 1.0 × 10−2 | 4.0 × 10−5 | ||||
Material constants
| Material | SS400 | Coil (copper) | EPE | Air |
|---|---|---|---|---|
| Electric conductivity σ (S/m) | 7.5 × 106 | No eddy current | ||
| Mass density ρ (kg/m3) | 7,850 | 8,880 | 1,800 | 1.18 |
| Heat capacity c (J/kg m3) | 465 | 386 | 2,302 | 1,007 |
| Thermal conductivity λ (W/m K) | 43 | 398 | 0.067 | 0.0261 |
| Viscosity η (Pa s) | Solid | 1.862 × 10−5 | ||
| Coefficient of volumetric expansion β (1/K) | 3.32 × 10−3 | |||
4.2 Coupled analysis with thermal fluid and heat conduction
The induction heating apparatus is analyzed by the coupled electromagnetic-field with heat conduction and thermal-fluid analysis. Figure 7(a) shows the result obtained by the experiment at time 180 s. Figure 7(b) shows the temperature distribution obtained by FE–FD–HS at time 180 s.
Comparison of the experiment and numerical simulation. (a) Experiment (t = 180 s), (b) numerical simulation (t = 180 s).
Figure 8 shows the advection effect of heat due to natural convection. A numerical evaluation of the temperature-rise characteristics at points α and γ shown in Figure 4 is shown in Figure 9.
Advection effect of heat by natural convection. (a) t = 45 s, (b) t = 90 s, (c) t = 135 s, (d) t = 180 s.
The numerical evaluation temperature characteristics of the heated target at points α and γ.
FE–FD–HS can follow experiment result, and the validity is confirmed by comparing the results of FE–FD–HS with those from the experiment. Because the discrepancy between FE–FD–HS and FE–FD at time 180 s is about 0.8 K at point α and γ, it is confirmed that heat-transfer effect of the heat target is affected by natural convection. The discrepancy between FE–FD–HS and FE–FD is slightly larger at point γ, because the fluid near point γ is cooler than the fluid near point α.
4.3 Heat-conduction analysis using local heat transfer coefficient
Figure 10 shows the distribution of the local heat transfer on the heated target. Its distribution is derived from the coupled analysis between the electromagnetic field (FEM), heat conduction (FDM), and thermal fluid (HSMAC). It can be seen that the distribution of the heat-transfer coefficient is locally quite different owing to the development of natural convection. Therefore, in the heat conduction analysis using the heat-transfer boundary, the actual phenomenon cannot be modeled by setting the heat-transfer coefficient to a constant value. Because the discrepancy between the surface temperature and the fluid temperature is low, the heat-transfer coefficient near the center of the heated target is evaluated as higher. A numerical evaluation of the temperature-rise characteristics at points α and γ shown in Figure 4 is shown in Figure 11.
Local heat-transfer distribution. (a) t = 45 s, (b) t = 90 s, (c) t = 135 s, (d) t = 180 s.
The numerical evaluation temperature characteristics of the heated target at points α and γ.
In this study, LHTC refers to the localized heat-transfer coefficient, and FE–FE (using LHTC) is shown as the temperature-rise characteristic of the heat-conduction analysis using the local heat-transfer coefficient obtained by the thermal-fluid analysis. In addition, FE–FD and FE–FE represent heat-conduction analysis without the heat-transfer effect. Because the temperature reduction by thermal fluid analysis (difference between FE–FD–HS and FE–FD) and temperature reduction by heat-conduction analysis using the local heat-transfer coefficient (difference between FE–FE (using LHTC) and FE–FE) are the same, it is observed that the local heat-transfer coefficient is calculated accurately.
5 Conclusion
The local heat-transfer coefficient is calculated by the convection distribution, which is simulated from the coupled electromagnetic field with heat conduction and thermal-fluid analysis. As a result, it can be seen that the heat-transfer coefficient is locally distributed because it is greatly influenced by the behavior of natural convection. Therefore, setting the heat-transfer coefficient to a constant value on the surface of the heated target is not suitable for modeling the actual phenomenon. Furthermore, it is verified that the heat-conduction analysis using the local heat-transfer coefficient can reproduce the temperature reduction of the thermal fluid analysis.
References
[1] Preis K, Biro O, Buchgraber G, Ticar I. Thermal-electromagnetic coupling in the finite-element simulation of power transformers. IEEE Trans Magn. 2006;42(4):999–1002.10.1109/TMAG.2006.871439Search in Google Scholar
[2] Hirono K, Hoshino R, Kamiya T, Wakao S, Okamoto Y, Jeon W. Design optimization of primary core in induction heating roll by combination of 2D level-set method and 3D coupled magnetic-thermal FEM. IEEE J Ind Appl. 2018;7(1):64–72.10.1541/ieejjia.7.64Search in Google Scholar
[3] Kitak P, Glotic A, Ticar I. Heat transfer coefficients determination of numerical model by using particle swarm optimization. IEEE Trans Magn. 2014;50(2):7023104.10.1109/TMAG.2013.2282409Search in Google Scholar
[4] Koshizuka S, Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nuclear Sci Eng. 1996;123:421–34.10.13182/NSE96-A24205Search in Google Scholar
[5] Muramatsu K, Takahashi N, Mimura T. Magneto-thermal-fluid analysis taking account of natural convection using semi-Lagrange coordinate system. IEEE Trans Magn. 1999;35(3):1670–3.10.1109/20.767335Search in Google Scholar
[6] Yamamoto T, Matsuzawa S, Ota T, Hirata K. Numerical analysis of ion behavior considering charging effect of a dielectric body. IEEE Trans Magn. 2017;53(6):7205004.10.1109/CEFC.2016.7816264Search in Google Scholar
[7] Jiangfei L, Long J, Lian Y, Zhizhong F, Bo L, Wenxue C. Comparison of finite difference and finite volume method for numerical simulation of driven cavity flow based on MAC. Int Conf Comp Inf Sci. 2013;891–4.10.1109/ICCIS.2013.239Search in Google Scholar
[8] Itu LM, Suciu C, Moldoveanu F, Postelnicu A. Optimized GPU based simulation of the incompressible Navier–Stokes equations on MAC grid. 10th RoEduNet Inter. Conf., Iasi, Romania, 5–8 (2011).10.1109/RoEduNet.2011.5993692Search in Google Scholar
[9] Hirt C-W, Nichols B-D, Romero N-C. SOLA-A numerical solution algorithm for transient fluid flows. Los Alamos Scientific Laboratory Report, LA-5852, Add, 1975.10.2172/4205348Search in Google Scholar
[10] Marignetti F, Colli VD, Coia Y. Design of axial flux PM synchronous machines through 3-D coupled electromagnetic thermal and fluid-dynamical finite-element analysis. IEEE Trans Magn. 2008;55(10):3591–601.10.1109/TIE.2008.2005017Search in Google Scholar
[11] Klomberg S, Farnleitner E, Kastner G, Biro O. Comparison of CFD analyzing strategies for hydrogenerators. Int Conf Electr Mach. 2014;1990–1995.10.1109/ICELMACH.2014.6960457Search in Google Scholar
[12] Fujiwara K, Nakata T, Ohashi H. Improvement of convergence characteristic of ICCG method for the A-ϕ method using edge elements. IEEE Trans Magn. 1996;32(3):804–7.10.1109/20.497363Search in Google Scholar
[13] Leonard BP. A survey of finite differences with upwinding for numerical modelling of the incompressible convective diffusion equation. Comput Techniq Trans Turbul Flow. 1981;2:1–35.Search in Google Scholar
© 2020 Satoshi Namiki et al., published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.
