CFD Analysis of Natural Oil Convection in Three Phase Transformers between Detailed and Simplified Winding Models ()
1. Introduction
Transformer plays a vital role in modern power systems, not only in construction but also in day-to-day operation. Their basic function is made possible by electromagnetic principles, where the windings are magnetically linked through a core [1]. The health and lifetime of a transformer largely depend on the thermal aging of its insulation system, particularly the oil-impregnated paper. As the transformer operates, winding losses generate heat, and if not managed properly, this accelerates insulation degradation. Typically, transformers rely on natural oil circulation between the windings and oil ducts, with heat eventually transferred to radiators through natural or forced air convection. Designing effective cooling systems, however, is not straightforward. Engineers must carefully decide the type and size of cooling equipment, such as radiators, to prevent overheating. Early in the design stage, fast yet reliable estimates are needed so that a transformer can be proposed within cost and time constraints. Conventional methods, such as thermal network models (THNM), provide average temperature predictions, but they fall short in showing detailed temperature variations and hot-spot locations. For instance, Wang et al. [2] present that conventional THNM can falsely identify the hottest disc temperature in winding assemblies compared to computational fluid dynamic (CFD). Thus, these limitations make it difficult to ensure safe and efficient designs.
To address these challenges, engineers are increasingly turning to CFD. CFD is a powerful simulation tool that allows detailed analysis of heat transfer, oil flow velocity, and temperature distribution within transformers [3]. Unlike traditional methods, CFD can provide a clear picture of where heat builds up, how oil circulates, and how effective the cooling system is. Although CFD has become more popular in recent years, it comes with challenges of its own. The simulations require specialized knowledge, significant computing resources, and time, which makes them less practical for routine design tasks [4]-[7].
Research in this area has demonstrated the effectiveness of CFD in transformer analysis. For example, three-dimensional (3D) CFD models have been used to predict temperature distribution in windings more accurately than traditional thermal models [8] [9]. Studies have also explored issues such as oil leakage, nanoparticle-enhanced fluids, and the impact of manufacturing imperfections using CFD [10]-[12]. These efforts highlight CFD’s role in identifying hot-spot temperatures, improving cooling design and consequently extending the transformer’s lifespan.
This paper aims to explore the use of CFD for both simplified winding model (SW) and detailed winding model (DW) in transformer design. Specifically, it compares oil flow velocity and thermal distribution patterns using two-dimensional (2D) CFD simulations. The study also examines whether SW can be useful during the early stage of design process, where not only fast but reliable results are needed. Two transformers are used as case studies: a 1000 kVA, 11.0/0.415 kV unit and a 30 MVA, 33/11 kV unit. The 1000 kVA transformer uses copper foil conductors with layer-type windings, while the 30 MVA transformer uses rectangular conductors with helical-disc windings. By analyzing these two cases, the study provides practical insights into how CFD can improve transformer cooling design and thermal performance prediction.
2. Modeling and Simulation Setup
2.1. Material Properties
For this study, both the oil and solid materials of the transformers were assigned specific thermal and physical properties. Since transformer oil plays a critical role in heat transfer, its properties were modeled as varying with temperature, while pressure was assumed to remain constant. The transformer oil temperature-dependent relations used in this simulation can be found in Equations (1) to (4) [2]:
(1)
(2)
(3)
(4)
where
is the oil density,
is the thermal conductivity of oil,
is the oil specific heat capacity at constant pressure,
is the oil dynamic viscosity, and
is the oil temperature. The solid materials were treated as isotropic, meaning their properties were considered uniform in all directions and independent of changes in temperature or pressure. For the copper conductor, the density was set to 8933 kg·(m3)−1, the specific heat capacity to 385 J·kg−1·K−1, and the thermal conductivity to 401 W·m−1·K−1. The insulation paper was assigned a density of 930 kg·(m3)−1, a specific heat capacity of 1340 J·kg−1·K−1, and a thermal conductivity of 0.19 W·m−1·K−1.
By defining these properties, the thermal performance of both oil and solid components can be represented realistically, ensuring that the simulation results closely reflect actual transformer behavior.
2.2. Small Distribution Transformer of 1000 kVA 11/0.415 kV
The transformer is considered a hermetically sealed distribution-type unit operating with a single cooling mode of ONAN (Oil Natural Air Natural). The low-voltage (LV) winding, rated at 0.415 kV, was constructed from 16 turns of copper foil measuring 430 mm × 1.3 mm, separated by 0.07 mm diamond-dotted insulation paper. To improve cooling, an axial oil duct of 4.5 mm pressboard was introduced after the eighth turn. There are two parts to the LV winding: S1 (inner coil) and S2 (outer coil). This design ensures that cooling oil naturally flows from the bottom of the tank, through the axial duct, and exits at the top of the winding. The winding arrangement is illustrated in Figure 1. When fully loaded, the LV winding carries a phase current of 1333.33 A, which corresponds to a current density of 2.38 A·mm−2 in the LV winding. The total winding losses are made up of resistive and eddy current losses. The resistive (ohmic) loss,
is given by Equation (5).
(5)
where J is the current density, σ is the conductivity, L is the length,
is the conductor density. The expressions for the axial
and radial
eddy losses can be found in Equations (6) and (7) respectively.
(6)
(7)
By combining resistive, axial and radial eddy losses, the total LV winding loss per phase is estimated at 1036 W.
Figure 1. Close-up of the 2D axis-symmetric winding in detailed 0.415 kV LV winding model.
2.3. Medium Power Transformer of 30 MVA 33/11 kV
This transformer is a conservator-type with ONAN/ONAF (Oil Natural Air Natural/Oil Natural Air Force) cooling modes. The 11 kV low-voltage (LV) winding is a type of helical winding comprising 95 number of turns, with each turn consisting of 24 parallel strands. To enhance thermal performance, axial cooling ducts are introduced after every 12 conductors using a cooling band of 5 mm. The horizontal cooling ducts are formed by 3 mm axial spacers positioned between discs. The winding model is sectionized into two parts: S1 (inner coil) and S2 (outer coil). Detailed winding configuration follows a bottom to top oil flow velocity pattern as depicted in Figure 2. Under rated operating conditions, the full load current of LV winding is 1575 A, corresponding to a current density of 2.22 mm−2. Using the established formulations Equation (5) to (7), the total winding losses are calculated to be approximately 16.36 kW per phase.
![]()
Figure 2. Close-up of the 2D axis-symmetric winding in detailed 11 kV LV winding model.
3. Methodology
This study aims to develop a reliable method for predicting thermal distribution in transformer windings using CFD. For this reason, a simplified 2D winding model has been constructed and investigated. As seen in Figure 3(b) and Figure 3(f), these winding models employed a single big conductor representing all individual conductors as shown in Figure 3(a) and Figure 3(e). This approach is justified by the significant difference in thermal conductivity between copper and Kraft paper. Since copper conducts heat far more efficiently than the insulation, the temperature gradient across the conductor is insignificant and can be disregarded [13].
Figure 3. Detailed model (DW) versus Simplified models (SW) of low voltage winding (a) 1000 kVA—DW (b) 1000 kVA—SW (c) 1000 kVA Mesh plot—DW (d) 1000 kVA Mesh plot—SW (e) 30 MVA—DW (f) 30 MVA—SW (g) 30 MVA Mesh plot—DW (h) 30 MVA Mesh plot—SW.
It is observed that the conductor cross-sectional area increases progressively as can be seen from Figure 3(a) to Figure 3(b) and from Figure 3(e) to Figure 3(f). To maintain consistent cooling performance across conductors, it is important that each disc has the same effective contact area as the surrounding oil. This relationship can be described using Newton’s Law of Cooling [13]:
(8)
Here, is the heat dissipated by the system, is the heat transfer coefficient, and represents the surface and oil temperatures, respectively. Term A refers to the cooling surface area of a single turn or disc in contact with oil, which can be expressed as:
(9)
In this equation, AD is the axial and RD is the radial dimensions of conductors, while OD and ID are measured for winding inner and outer diameters respectively. The average circumferential length is denoted by lmean and determined by:
(10)
V is the winding volume, can be computed using:
(11)
Here, RR refers to the radial dimension of the winding. The overall heat generation, Qs within the winding is represented in terms of volumetric heat loss such as:
(12)
where Ptot is the total winding loss. This stepwise formulation ensures that the influence of conductor geometry and oil contact area on the thermal behavior of the winding is fully captured.
According to Santisteban et al. [13], the heat transfer formulation in Equation (13) is employed to ascertain the temperature distribution within the windings and the surrounding oil, while Equation (14) is the Navier-Stokes and Equation (15) is the conservation of mass for fluid flow equations are utilized to analyze the oil flow velocity problem respectively. The expression for these governing equations is:
(13)
(14)
(15)
where u is the vector of velocity field, p is the pressure and g stands for the gravitational acceleration.
Variations in buoyancy forces brought on by changes in oil density because of thermal expansion are considered by the buoyancy component on the right side in Equation (14). Under steady-state conditions, the first term of transient component in the Equation (13) can be neglected. Additionally, the heat generation term does not apply to solid materials such as copper windings and insulation paper.
Solerno et al. [14] explains the role of three dimensionless factors which are commonly used to describe the type of oil flow velocity within the windings. Firstly, the Reynolds number (Re) that characterizes whether the flow is laminar, transitional or turbulent. The flow remains laminar if the Re is below 2000. Second is the Grashof number (Gr), which reflects the role of natural convection by comparing buoyancy to viscous forces. Thirdly, the Prandtl number (Pr) describes the relative efficiency of momentum and energy transfers by diffusion in the fluid flow and thermal boundary. The ratio Gr/Re2 provides insights into the dominant flow regime. Specifically, Gr/Re2 1 indicates natural convection, while Gr/Re2 1 corresponds to forced convection. These dimensionless factors are expressed as the following.
(16)
(17)
(18)
where U represents the fluid flow velocity, L is the length, β is the thermal expansion coefficient, ν denotes the kinematic viscosity and α is the fluid thermal diffusivity.
In the simulations, oil flow was assumed to remain laminar due to the fluid’s high viscosity. A two-dimensional axisymmetric winding model was adopted to reduce computational demand. Cylindrical pressboards were treated as adiabatic, justified by their low thermal conductivity (approximately 0.1 W·m−1·K−1) [15]. Uneven insulation surfaces were neglected with no-slip boundary condition was applied [6]. Although in practice eddy current losses vary across individual turns or discs, for simplicity, winding losses were uniformly distributed in the simulations.
Insulation paper and copper windings were modeled as homogeneous and temperature independent. During ONAN operation, oil flow was driven by buoyancy forces, with ambient condition kept constant and open boundaries applied at the inlet and outlet [6] [16].
4. Results and Discussions
4.1. 1000 kVA, 11/0.415 kV Small Distribution Transformer
Table 1 presents the simulation outcomes for a 1000 kVA transformer. The results show slight differences between DW and SW models. The average winding temperature (Tave) of the SW model was approximately 1.42% higher than that of the DW model, while the difference in the average winding temperature gradient (ΔTT-B) was just 0.55 K. Furthermore, the hot-spot temperature (Thot-spot) of the DW model was 1.46˚C lower compared to the SW model, and the maximum temperature difference (Tmax – Tmin) varied by merely 2.9%.
Table 1. Temperature data between DW and SW models.
|
Tave (˚C) |
ΔTT-B (˚C) |
Thot-spot (˚C) |
Tmax – Tmin (˚C) |
DW |
84.70 |
(89.94 − 76.92) = 13.02 |
90.44 |
(90.44 − 75.66) = 14.78 |
SW |
85.90 |
(91.17 − 77.60) = 13.57 |
91.90 |
(91.90 − 77.55) = 14.35 |
Figure 4 and Figure 5 illustrate the surface velocity and temperature distribution profiles for both models. The hot-spot region for DW and SW was located between 420 mm and 430 mm of winding height.
Figure 6 compares the mean oil velocity across vertical oil ducts in the two models. The oil flow velocity was measured horizontally at intervals of 50 mm across the winding height. Results indicate that at oil duct no.2, i.e. DW_C2 and SW_C2 exhibited the highest mean oil flow velocity for both DW and SW models. In contrast, the oil flow velocity in oil duct no. 1 and no. 3, i.e. DW_C1, DW_C3, SW_C1 and SW_C3 showed comparable values and similar flow patterns.
Figure 4. The magnitude of oil flow velocity (mm·s−1) for (a) DW, (b) SW.
Figure 5. Temperature distribution (˚C) plots on (a) DW, (b) SW.
Figure 6. Oil flow velocity inside cooling ducts between DW and SW with respect to Figure 4.
Figure 7. Winding temperature for every section between DW and SW models.
Figure 7 shows the average winding temperatures, which display only minor variation between the two models. The temperature distribution remained relatively uniform across the winding height. Notably, the temperature increased almost linearly with winding height before stabilizing toward the top of winding, approximately between 350 mm and 430 mm.
Figure 8. (a) DW-temperature distributions; (b) SW-temperature distributions; (c) DW-oil flow velocity magnitudes; (d) SW-oil flow velocity magnitudes.
Figure 8 presents the 3D plot comparisons between DW and SW, offering a clearer visualization of the flow and temperature distribution trends. These results complement the quantitative analysis summarized in Table 1, providing both numerical and graphical insights into model performance.
4.2. 30 MVA, 33/11 kV Medium Power Transformer
Table 2 summarizes the simulation results for a 30 MVA transformer with 94 discs. The SW model exhibited a winding temperature gradient (ΔTT-B) nearly 280% higher than DW, whereas its maximum temperature difference (Tmax – Tmin) was 14.49˚C lower. The hot-spot temperature of SW was also reduced by 12.89˚C compared to DW, largely due to the absence of insulation paper which has lowered thermal resistance. Despite these differences, average winding temperatures (Tave) varied only slightly between models.
Table 2. Temperature data between DW and SW.
|
Tave (˚C) |
ΔTT-B (˚C) |
Thot-spot (˚C) |
Tmax – Tmin (˚C) |
DW |
66.82 |
(64.72 − 62.69) = 2.03 |
84.20 |
(81.99 − 56.10) = 25.89 |
SW |
64.84 |
(35.59 − 57.88) = 7.71 |
65.73 |
(69.10 − 57.66) = 11.44 |
Figure 9 shows the oil velocity and temperature distribution. Both models recorded a similar maximum oil velocity of about 51 mm·s−1. Hot spots occurred on disc no. 93 in DW at 84.82˚C and disc no. 94 in SW at 65.73˚C. In SW disc no.49 to 51 showed higher temperature which exceeded 80˚C. For DW, lower pressure differences at disc junctions caused regions of weak oil circulation within horizontal ducts.
Figure 10 and Figure 11 present a comparison of the models with respect to oil flow velocity and thermal distribution, respectively. The oil flow velocity was determined as the mean flow magnitude measured at the midpoint of each horizontal oil duct, whereas the temperature was expressed as the mean disc temperature.
As seen in Figure 10, the oil flow velocity profiles of SW_S1 and SW_S2 exhibit similar overall behaviour, with consistent variations across most oil ducts, except at duct no.1 and duct no.24. Despite this similarity, SW_S1 recorded the highest velocity of 1.5 mm·s−1 at duct no.1, while SW_S2 reached its peak at 0.76 mm·s−1 in duct no.33.
For the DW model, the flow distribution was nearly identical for S1 and S2. In both cases, the velocity gradually increased from duct no. 19 to duct no.27 before decreasing sharply to nearly zero at duct no.50 (S2) and duct no.53 (S1). Beyond this point, the flow rose steeply, peaking at 2.43 mm·s−1 and 2.49 mm·s−1 for S1 and S2, respectively. At duct no.93, the velocity again decreased, settling at 0.27 mm·s−1 for S1 and 1.09 mm·s−1 for S2.
From Figure 11, both models demonstrate generally comparable thermal distribution patterns across the discs. In the SW model, the temperature rises in an almost linear trend for both S1 and S2. At disc no.1, S1 measured at 53.72˚C rises to 64.03˚C at disc no. 94, likewise for S2 but with slightly higher temperature of 65.73˚C at disc no.94. In contrast, the DW model shows a less uniform distribution, characterized by three distinct hot-spot regions located at discs no.10 (S2), 50 (S2), and 93 (S1). The corresponding mean temperatures at these positions were 71.62˚C, 76.54˚C, and 78.98˚C respectively. The maximum temperature of 84.24˚C arose at disc no.93 (S1), while the minimum value of 57.98˚C was observed at disc no.1 (S2).
Figure 9. Two-dimensional plots: (a) DW—Magnitude of oil flow velocity (mm·s−1) (b) SW—Magnitude of oil flow velocity (mm·s−1) (c) DW—thermal distribution (˚C) (d) SW—thermal distribution in (˚C).
Figure 10. Oil flow velocity inside 3 mm height of radial cooling ducts.
Figure 11. Average winding temperature comparison between SW and DW.
Figure 12 highlights the hot-spot region in the DW model at disc no.93. This location coincides with the lowest resultant radial oil flow velocity, as indicated by the black circle. As illustrated in Figure 13(b), the radial component of oil flow velocity (ranging between +0.6 mm·s−1 and −0.6 mm·s−1) move at nearly equal magnitudes but in opposite directions. This phenomenon produces oil stagnation or weak circulation zones which are similar to eddies.
Figure 12. Close-up: (a) hot-spot temperature location (b) oil flow velocity in radial component (mm·s−1).
Consequently, heat is not effectively transported away from this region, leading to the formation of the hottest zone. The observed flow distribution of the oil therefore provides a direct explanation for the occurrence and positioning of this winding hot-spot.
Figure 13 displays the 3D plots for the oil velocity and thermal distributions of SW and DW, respectively. The analysis of Figure 13 was tabulated in Table 2.
Figure 13. (a) Distribution of temperature for DW, (b) Distribution of temperature of SW, (c) Oil flow velocity of DW and (d) Oil flow velocity of SW.
5. Conclusion
This study showed that buoyancy-driven oil flow velocity is strongly dependent on gravitational and temperature effects, making flow direction and hot spot location difficult to predict. For 1000 kV distribution transformers with copper foil windings, SW and DW models produced similar thermal results. In contrast with 30 MVA transformers with disc winding, the SW approach failed to capture DW behavior in both thermal distribution and oil flow. Therefore, simplification is acceptable for small distribution transformers but not for medium and large power transformers, where detailed winding models are essential.
Acknowledgements
The authors wish to acknowledge the Universiti Teknikal Malaysia Melaka (UTeM) and Centre for Research and Innovation Management (CRIM) for funding this research paper.