CHEMICAL ENGINEERING TRANSACTIONS VOL. 81, 2020 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Petar S. Varbanov, Qiuwang Wang, Min Zeng, Panos Seferlis, Ting Ma, Jiří J. Klemeš Copyright © 2020, AIDIC Servizi S.r.l. ISBN 978-88-95608-79-2; ISSN 2283-9216 Multi-Physics Coupling Simulation of Heat Transfer in Transformer Winding Xiaoyu Jiaa, Wenrong Sic, Chenzhao Fuc, Yining Wub, Qiuwang Wangb, Mei Lina, Jian Yangb,* aSchool of Energy and Power Engineering, Xi'an Jiaotong University, Xi'an 710049, China bKey Laboratory of Thermo-Fluid Science and Engineering, MOE, School of Energy and Power Engineering, Xi'an Jiaotong University, Xi'an, Shaanxi, 710049, P.R. China cState Grid Shanghai Electrical Power Research Institute, Shanghai 200437, PR China yangjian81@mail.edu.cn In the electric power industry, the requirements for capacity, voltage level and compactness of transformers are getting higher and higher, which lead to significant problems, such as increased loss and more serious hot spots phenomena. Based on the three-dimensional finite element method, this paper establishes the bidirectional coupling between electromagnetic field, flow field and temperature field to calculate the temperature rise of 400 kVA oil-immersed transformer at full load and compare it with the results obtained by the heat-flow unidirectional coupling method (constant loss). The results show, due to the interaction between electromagnetic field and temperature field, the coil loss density increases, but the core loss density does not change. The hot spot temperature of coil increases and the location of coil hot spot changes. The temperature and position of the core hot spot do not change. The loss density of coil decreases as the inlet velocity increases, the hot spot temperature of coil drops faster. When the inlet velocity exceeds 0.5 m/s, increasing the inlet velocity can no longer effectively promote heat exchange. With the increase of the inlet velocity, the interaction between electromagnetic field and temperature field is weakened, but it still has a great impact on the loss density of coil. The multi-physics coupling method is of great significance to improve the accuracy of transformer temperature rise calculation. 1. Introduction Power transformers are one of the most critical devices in power system (Zhou et al., 2017), its insulation performance is related to the heat generation and heat dissipation of the winding. If the winding operating temperature is too high, it will cause serious electrical or mechanical failure (Zhang, 2017). Accurately calculating the temperature rise and predicting the temperature and location of hotspot is of great significance to improve the safety and stability of transformer operation and extend its service life (Acharya et al., 2017). Deng et al. (2019) proposed a prediction method to predict the hot spot temperature of a 10 kV oil-immersed transformer, took feature points and transformer load rate as the input parameters of a machine learning model established by support vector regression (SVR) to describe their relationships. Gamil et al. (2018) established a theoretical thermal model based on thermal hydraulic principles to calculate the temperature of winding and inserted the data from cooling-curves of a set of different transformer designs to correct the calculated gradient and recalculate the oil temperature at top. Numerical simulation methods to analyze electrical equipment temperature distribution have been widely used. Sitar et al. (2017) calculated the tank loss by the edge element version of T-Ω method, then used it as a constant heat source to conduct a temperature FEM calculation. Liu et al. (2019) established a two-dimensional model of oil-immersed transformer, used the finite element method to obtain the transformer winding loss, and then used it as a constant heat source to calculate the oil and temperature distribution of the transformer. Yang et al. (2019) adopted the principle of constant joule heating to transform the fluctuating current into a steady-state DOI: 10.3303/CET2081057 Paper Received: 06/04/2020; Revised: 18/05/2020; Accepted: 22/05/2020 Please cite this article as: Jia X., Si W., Fu C., Wu Y., Wang Q., Lin M., Yang J., 2020, Multi-Physics Coupling Simulation of Heat Transfer in Transformer Winding, Chemical Engineering Transactions, 81, 337-342 DOI:10.3303/CET2081057 337 current to obtain the internal dielectric loss of transformer bushing, then used it as a constant heat source to calculate the transient temperature field of bushing. Temperature rise of transformer is a multi-physics coupling process, but judging from the studies mentioned above, scholars currently calculate transformer temperature rise based on the assumption of constant losses, using the heat-flow unidirectional coupling method, which fails to consider the effect of temperature on losses. And scholars usually use the two-dimensional long and short axis slice model of the winding to replace the three-dimensional model, which cannot accurately predict the hot spot location of winding. This paper establishes the bidirectional coupling between electromagnetic-heat-flow physical fields for the first time to calculate the winding loss and temperature distribution of the transformer at full load and compares it with the results when the winding loss is constant to analyze the interaction between multi-physics. This work will contribute to the development of numerical calculation of transformers. 2. Physical model and calculation method 2.1 Physical model and geometric parameters 400 kVA-15 kV/400 V oil-immersed transformer is selected as the calculation model. The cooling method is forced oil circulation. It uses an oil pump. The power is greater than the natural oil circulation method, the cooling effect is better. The heat exchange effect can be adjusted by the inlet oil flow rate. Its physical model and size are shown in Figure 1. Including: transformer oil tank (Figure 1a) and winding (Figure 1b). To simplify the model, ignore the transformer oil tank thickness. Inlet velocity is 0.1 m/s. (a) (b) Figure 1: Transformer model and geometric parameters (a) Transformer tank and (b) Winding (core and coil) The physical properties of winding and transformer oil vary with temperature. Transformer oil physical properties are given by COMSOL finite element analysis software. The thermal conductivity and specific heat capacity of winding are show in Eq(1) - (4) (Tsili et al., 2012). Where cp, iron (T) and cp, copper (T) are the constant pressure heat capacity of core and coil, (W/(m·K)); λiron (T) and λcopper (T) are the thermal conductivity of core and coil, (J/(kg·K); T is temperature, K. ( ) ( ) ( ) 2-5 iron 8.64 10 - 273.15 - 0.104 - 273.15 404.18T T T =    + (1) ( ) ( ) ( ) 2-4 p, iron 2.91 10 273.15 0.522 273.15 431.88c T T T= −   − +  − + (2) ( ) ( ) ( ) 2-4 copper 1.22 10 - 273.15 - 0.128 - 273.15 83.71T T T =    + (3) ( ) ( ) 2-4 p, copper -3.20 10 - 273.15 0.221 ( - 273.15) 376.98c T T T=   +  + (4) σ (T) is the electrical conductivity of coil, S/m; μr is the relative permeability of core, as show in Eq(5)-(6). ( ) ( )( )0 ref1 / 1 -T T T  = + (5) 950 30 r = -j j    = − (6) The electrical conductivity of coil changes with temperature, the reference resistivity ρ0 is 1.6679 × 10-8 Ω×m, the resistivity temperature coefficient α is 3.9 × 10-3 1/K, the reference temperature Tref is 293.15 K. The relative 338 permeability of core is expressed in the plural form, μ is the real part, represents the magnetic storage capacity, μ is the imaginary part, represents the magnetic loss capacity, j is the imaginary unit of plural. 2.2 Control equations and calculation methods The transformer model is divided into a solid domain (winding) and a fluid domain (transformer oil). The electromagnetic-heat coupling equations of the solid domain are show in Eq(7)-(9). The heat transfer equation of winding Eq(7) and the electromagnetic filed equation Eq(9) are coupled by the magnetic loss equation Eq(8). The tank wall is magnetically insulated. ( )( ) 0eT T Q  + = (7) ( ) ( )* * 1 1 2 2 e rh ml Q Q Q Re J E Re j B H= + =  +  (8) ( ) ( )2 1 0A j T A A   −− + +  = (9) The constitutive equations of electromagnetic field are show in Eq(10)-(13), A B = (10) E j A= − (11) B H= (12) ( )J T E= (13) The coupled heat-flow equations of the fluid domain are show in Eqs(14) - (16), the mass conservation equation Eq(14), the oil heat transfer equation Eq(15) and the momentum conservation equation Eq(16) are coupled by the thermal properties of transformer oil (oil properties are given by COMSOL finite element analysis software). 0v = (14) ( ) ( ) ( ) ( ) [( ) ( )]t p T T T v T T c T      =   +   (15) ( ) ( )( ) [( )( ( ) )]TtT v v p T v v   = − + +  +  (16) Where λ (T) is the thermal conductivity, W/(m2·K); Qe is the winding loss density, Qrh and Qml are the electrical loss density and magnetic loss density, W/m3; Re is the real part of the imaginary number; J is the current density, A/m2; * E and * H are conjugate complex numbers of E and H ; B is the magnetic induction, T; A is the magnetic vector potential, Wb/m; E is the electric field strength, N/C; H is the magnetic field strength, A/m; j is an imaginary unit; ω = 2πf, is the phase angle; f is the electric excitation frequency, 50 Hz; γ is the dielectric constant, F/M; μ is the permeability, H/m; σ (T) is the electrical conductivity, S/m; v is the speed vector of transformer oil, m/s; ρ (T) is the transformer oil density, kg/m3; cp (T) is the constant pressure specific heat capacity of transformer oil, J/(kg·K); μt is the turbulent viscosity of transformer oil, kg/(m·s); σT is the Prandtl number of the transformer oil heat transfer equation; p is pressure, Pa; μ (T) is the dynamic viscosity of transformer oil, kg/(m·s). The boundary condition of the contact surface between the winding and transformer oil is fluid-solid coupling heat transfer boundary. The heat exchange between the tank wall and the outside air is natural convection, and the ambient temperature, Tamb is 293.15 K. The heat transfer at the top and bottom tank wall can be considered as natural convection heat transfer outside a horizontal wall, the heat transfer at the side wall can be regarded as the natural convection heat transfer outside a vertical wall. The governing equations are solved with the finite element analysis software COMSOL Multiphysics 5.2. The iterative solver GMRES of frequency domain-steady state solution is used to calculate the electromagnetic loss, temperature and oil distribution of the transformer at full load (COMSOL Multiphysics 5.2, 2016). When the calculated residual is less than 10-3, the calculation is considered to have reached convergence. 339 3. Grid independence test and model validations The memory resources required for multi-physics direct coupling method are extremely large, too fine meshing makes calculation impossible. Based on Richardson's extrapolation method (Roache, 1994), the grid convergence index (GCI) is used to verify the independence of the grid. The total number of three sets of grids is 395,980, 688,181, and 1,109,682. The GCI criterion number is 0.4 %. It can be considered that a grid- independent solution is obtained when the number of grids is 1,109,682. At present, no scholars have used the electromagnetic-heat-flow bidirectional coupling method to calculate the temperature rise of transformer. This paper checks each physics calculation model separately. When the electrical conductivity is constant (σ = 6 × 107 S/m), open and short circuit tests are performed on the transformer to check the electromagnetic loss model. The core loss obtained by the simulated open circuit test is 1,577.3 W, while the core loss obtained by Steinmetz formula (Kulkarni and Khaparde, 2017) is 1,601.74 W. The coil loss obtained by the simulated short circuit test is 2,269.93 W, and the DC loss of the coil according to Joule's law is 2,277.73 W. The results show that the electromagnetic loss calculation model is reliable. When studying the heat transfer characteristics of pipe cable in the previous article (Fu et al., 2018), a unidirectional coupling model of electromagnetic-heat-flow had been established. Following the previous study, an electromagnetic- heat-flow bidirectional coupling is established based on the change of electrical conductivity with temperature in the present paper. 4. Results and discussion When using the multi-physics coupling method, the electrical conductivity is affected by temperature, the high voltage coil loss density is 28,343 W/m3, the low voltage winding loss density is 39,516 W/m3, the core loss density is 10,697 W/m3. The temperature and oil distribution obtained by the multi-physics coupling method is shown in Figure 2. It shows that, there is a vortex at the bottom of the tank, and the oil flow velocity on both sides of the iron core is larger (Figure 2c). The core temperature distribution decreases from top to bottom (Figure 2b). The coil hot spot temperature is 371.06 K, which is located at 82.5 % of the total height of A-phase high voltage coil. The core hot spot temperature is 337.53 K, which is located on the upper yoke of core (93 % of the total height). (a) (b) (c) Figure 2: Temperature distribution of (a) Coil, (b) Core and (c) Speed distribution of oil in z direction (x-z section) obtained from the multi-physics coupling method When the electrical conductivity does not change with temperature (σ = 6 × 107 S/m), the winding loss is constant, the loss density of high and low voltage coils are 21,849 W/m3 and 31,181 W/m3, and the core loss density is 10,752 W/m3. The temperature and oil distribution obtained by the heat-flow unidirectional coupling method is shown in Figure 3. The coil hot spot temperature is 352.28 K, which is located at 87.5 % of the total height of B- phase low voltage coil. The core hot spot temperature is 337.66 K, which is located at the upper yoke of core (95.2 % of the total height). Compare the results of two methods, when the multi-physics coupling method is used, due to the electrical conductivity of coil changes with temperature, the loss density of high and low voltage coils increase by 29.72 % and 26.73 %. The coil hot spot temperature increases by 18.78 K, the hot spot location also changed. The core loss density is related to the magnetic flux density in the core, when the transformer is fully loaded, the magnetic flux density reaches saturation, and the core loss density and hot spot are basically unchanged. 340 (a) (b) (c) Figure 3: Temperature distribution of (a) Coil, (b) Core and (c) Speed distribution of oil in z direction (x-z section) obtained from the heat-flow unidirectional coupling method (σ=6×107 S/m, constant heat source) At different inlet velocity, the coil loss density obtained by the multi-physics coupling method and the percentage increase compared to the constant loss (σ = 6 × 107 S/m, the loss density of high and low voltage coils are 21,849 W/m3 and 31,181 W/m3) are shown in Table 1. With the increase of the inlet velocity, the coil loss density and the percentage increase of coil loss density decreases. Table 1: Coil loss density obtained by the multi-physics coupling method and percentage increase compared to constant loss Inlet velocity (m/s) Loss density of High voltage coil (W/m3) Loss density of low voltage coil (W/m3) Percentage increase of high voltage coil Percentage increase of low voltage coil 0.1 28,343 39,516 29.7 % 26.7 % 0.2 25,745 36,021 17.8 % 15.5 % 0.3 24,795 34,730 13.5 % 11.4 % 0.4 24,280 34,031 11.1 % 9.1 % 0.5 23,948 33,588 9.6 % 7.7 % 0.6 23,716 33,278 8.5 % 6.7 % The variation of hot spot temperature with inlet velocity is shown in Figure 4. As the inlet velocity increases, the loss density of coil decreases, the hot spot temperature of coil drops faster. The core loss density is unchanged when the transformer is fully loaded, so the hot spot temperature of core at different inlet velocity obtained by two methods is the same. When the inlet velocity exceeds 0.5 m/s, increasing the inlet velocity unable to further promote the disturbance of oil flow, and the temperature difference between winding and oil is small, so it can no longer effectively enhance heat transfer. 0.1 0.2 0.3 0.4 0.5 0.6 300 310 320 330 340 350 360 370 380 H o t s p o t te m p e ra tu re , T /( K ) Inlet velocity, V/(m/s) coil (multiphysics coupling) core (multiphysics coupling) coil (constant heat source) core (constant heat source) Figure 4: Variations of hot spot temperature with inlet velocity As can be seen from Table 1 and Figure 4, as the inlet velocity increases, the interaction between electromagnetic field and temperature field weakens, but the electromagnetic loss is still significantly affected by temperature. When the heat dissipation of the winding is poor or the transformer is a large-capacity transformer with large loss, it will also have a greater impact on the hot spot temperature. 341 5. Conclusions Compare the results of electromagnetic-heat-flow bidirectional coupling method with the results of heat-flow unidirectional coupling method, the main findings are as follows: (1) When the inlet velocity is 0.1 m/s, the loss density of high and low voltage coils increase by 29.72 % and 26.73 %, the hot spot temperature of coil increases by 17.5 K, the hot spot location of coil also changed. (2) The loss density of coil decreases with the increase of inlet velocity, the hot spot temperature of coil drops faster. When the inlet velocity exceeds 0.5 m/s, increasing the inlet velocity can no longer effectively enhance heat transfer. (3) As the inlet velocity increases, the interaction between electromagnetic field and temperature field is weakened, but it still has a significant impact on the electromagnetic loss, and when the heat dissipation of the winding is poor or the transformer capacity is large, it will also affect hot spot. This article is of great significance to analyze the physical mechanism of multi-physics strength coupling and quickly predict the hot spot of transformers, but due to the simplification of the physical model, the results in this article may deviate from the actual situation. In the future, the porous media model will be applied to the coil, the stray losses on the metal components will be considered. Acknowledgement This work was supported by the S&T of State Grid Corporation of China Research and application of key technologies of integrated ultrasonic and UHF PD coupling and sensing in transformer under grant of 52094019000S. References Acharya S., Tapre P.C., 2017, Life assessment of transformer using thermal models, 2017 International Conference on Energy, Communication, Data Analytics and Soft Computing (ICECDS 1-2 Aug, 2017), Chennai, India, 3515-3520. COMSOL Multiphysics® v. 5.2., 2016. COMSOL AB, Stockholm, Sweden, , accessed 20/05/2020. Deng Y.Q., Ruan J.J., Quan Y., Gong R.H., Huang D.D., Duan C.H., Xie Y.M., 2019, A method for hot spot temperature prediction of a 10 kV oil-immersed transformer, IEEE Access, 7, 107380-107388. Fu C.Z., Si W.R., Quan L., Yang J., 2018, Numerical study of convection and radiation heat transfer in pipe cable, Mathematical Problems in Engineering, 11, 1-12. Gamil A., Al-Abadi A., Schatzl F., Schluecker E., 2018, Theoretical and empirical-based thermal modelling of power transformers, IEEE International Conference on High Voltage Engineering and Application (ICHVE 10-13 September 2018), Athens, Greece, 1-4. Kulkarni S.V., Khaparde S.A., 2017, Transformer engineering: design, technology, and diagnostics (second edition), China Machine Press, Beijing, China. Liu G., Jing Y.J., Ma Y.Q., Sun L.P., Chi C., 2019, 2-D transient temperature field simulation of oil-immersed transformer based on hybrid method, High Voltage Apparatus, 55, 82-89. Roache P.J., 1994, Perspective: A method for uniform reporting of grid refinement studies, Journal of Fluids Engineering, 116, 405-413. Sitar R., Sulc I., Janic Z., 2017, Prediction of local temperature rise in power transformer tank by FEM, 4th International Colloquium on Transformer Research and Asset Management (10-12 May 2017), Pula, Croatia, 202, 231-239. Tsili M.A., Amoiralis E.I., Kladas A.G., Souflaris A.T., 2012, Power transformer thermal analysis by using an advanced coupled 3D heat transfer and fluid flow FEM model, International Journal of Thermal Sciences, 53, 188-201. Yang Z.F., Ruan J.J., Huang D.C., Du Z.Y., Tang L.Z., Zhou T.T., 2019, Calculation of hot spot temperature of transformer bushing considering current Fluctuation, IEEE Access, 7, 120441-120448. Zhou G.H., Liu G., Li M., 2017, Fault diagnosis of power transformer based on gas characteristic analysis, Chemical Engineering Transactions, 59, 793-798. Zhang K., 2017, Design of transformer oil temperature intelligent monitoring and alarm system based on PLC, Chemical Engineering Transactions, 62, 877-882. 342