Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 31, N o 3, September 2018, pp. 411-423 https://doi.org/10.2298/FUEE1803411S FEM CFD ANALYSIS OF AIR FLOW IN KIOSK SUBSTATION WITH THE OIL IMMERSED DISTRIBUTION TRANSFORMER  Stevan Stanišić 1 , Milica Jevtić 1 , Bhaba Das 2 , Zoran Radaković 1 1 Faculty of Electrical Engineering, University of Belgrade, Belgrade, Serbia 2 Engineering Department, ETEL Ltd, Auckland 0640, New Zealand Abstract. In practice of loading of oil-immersed distribution transformers, there is a need to have lumped thermal model, requiring no big computational resources and computational time. One such model is presented in international transformer loading guide (IEC 60076-7), where heat transfer inside the transformer is modeled. In case of indoor transformer operation, this model does not consider transient thermal phenomena in the room. We developed a lumped model that includes heat transfer in the transformer room. In scope of the research, we also built FEM CFD (finite element method, computational fluid dynamics) model of air flow and heat transfer. The purpose of FEM CFD was to make a better insight into air flow, i.e. to study the simplifications introduced in lumped model and suggest potential improvements. This paper presents results achieved with FEM CFD. The considered case was the transformer with natural oil and natural air flow (ONAN). Key words: Indoor transformer station, Thermal model, Finite Element Method, Computational Fluid Dynamics 1. INTRODUCTION It is well-known that the temperatures at the hottest position (hot-spot) in solid insulation and the hottest oil are the main factors which define possible transformer load (current) at specific ambient conditions (ambient temperature). The majority of the published research relate to the modeling of the heat transfer inside the transformer tank and from the tank and coolers (radiators) to the outer cooling medium. At this point, the losses depend on the load and the average winding temperature. So, the coupled calculation of the temperatures and the losses has to be performed. Nowadays, there is a strong demand for saving the space above the ground surface when placing the distribution substations in urban areas. On the other hand, establishing of air flow cooling the transformer can be more difficult and very restricted for underground placement. There is a need to have a calculation method for quantifying this effect. A common solution is prefabricated concrete transformer substation. So far, there Received September 29, 2017; received in revised form December 25, 2017 Corresponding author: Zoran Radaković Faculty of Electrical Engineering, University of Belgrade, 73 Kralj Aleksandar Blvd, 11000Belgrade, Serbia (E-mail: radakovic@etf.rs) 412 S. STANIŠIĆ, M. JEVTIĆ, B. DAS, Z. RADAKOVIĆ are compact substation designs. Technical issues for such substations include sizing of the ventilation openings and choosing the opening types, taking protection into consideration as well (safety of human beings and animals against touching the metallic parts under voltage, entry of the animals, rain etc.). This research was initiated with the task to optimize size and placement of ventilation openings for compact kiosk substation produced by ETEL Ltd, New Zealand. Ventilation opening sizing is an old engineering problem - some 40 years ago there were investigations about it [1, 2, 3], but in modern engineering practice (especially for smart grid concept) there is a need not only to perform the calculations in steady states, but also to develop the dynamic thermal models for the estimation of possible overloading. Actual loading guide IEC 60076-7 considers additional heating of the transformer due to its positioning in enclosure in simplified manner. The rated top oil temperature rise is corrected (increased) for the difference of air temperature in the enclosure minus ambient air temperature. The standard contains the typical values of this temperature increase for different types of the enclosure, the transformer rated power and the number of transformers in the enclosure. In our previous publication [4], we presented the dynamic thermal model for prefabricated concrete enclosure, typically used in Power Distribution Company, Belgrade, Serbia. That model was based on the empirical data. In our recent publication [5] we published more physical factor based model. As far as we know, no other research efforts similar to those presented in [5] have been made by other authors. Meanwhile, we improved the lumped model announced in [5] and will publish such an upgraded model, with additional on-site tests, in a future paper. The lumped model is suitable for on-line applications, but is simplified and of limited accuracy. In recent years, FEM CFD tools are becoming more and more present in the process of research and development of optimal cooling systems for power transformers. So far, main targets of FEM CFD analyses in this area are oil cooled core windings as well as both natural and forced transformer substation cooling. It is reported in [6] how finite element approach is applied on oil filled disc type windings for the purpose of locating hot spot. Cases regarding optimal design of indoor substation and ventilation are treated in [7, 8, 9, 10, 11]. There is even an example where FEM CFD tools were used to gain further insight in conjugate heat transfer for oil inside and air outside the radiators for both ONAN and ONAF transformers [12, 13]. This paper presents the results of the application of FEM CFD simulations, which gives detailed space distribution of air velocity and temperature. Thus, the simplifications in lumped model can be checked and lumped model potentially improved. The paper presents the experience about the application of FEM CFD simulations and their results. 2. GEOMETRY OF SIMULATED SUBSTATION The kiosk substation consists of transformer, HV and LV compartments. Figure 1 presents the geometry of the transformer compartment used in the model (base Sh = A x B = 1.56 x 1.32 = 2.06 m, height H = 1.33 m). The data about real kiosk transformer compartment ventilation openings is specified in Table 1. Table 2 contains the data about the simplified openings modeled as shown in Figure 1. Boundary condition with pressure head loss coefficient is assigned to the openings. The rated power of the transformer is 500 kVA and it is cooled by natural ventilation over the fins. FEM CFD Versus Lumped Thermal Model of Kiosk Substation with the Oil Immersed Distribution Transformer 413 Kiosk floor is 125 mm thick, kiosk walls containing the openings and kiosk ceiling are all 25 mm thick. In order to simplify the model, only kiosk transformer compartment interior surfaces were modeled, i.e. the wall thickness and the resistance to the heat conduction are neglected. Only the convection heat transfer is considered in the model. Convection heat transfer coefficient on the inner surfaces of the kiosk is determined from CFD calculation. Constant values of convection heat transfer coefficients on outer surfaces are determined using the equation from the theory of natural air flow near horizontal and vertical walls. The drawing of the modeled transformer compartment is shown in Figure 1, with openings protruding 25 mm to the outside, thus accounting all hydraulic resistances to air flow. Fig. 1 The geometry modeled with CFD FEM Table 1 Data about the real openings Openings Outlet Inlet Louvre Access Door (A) Louvre Side Panel (B) Cutout (D) Cutout (C) Hole Rows 10 10 7 18 Hole Columns 5 8 9 9 Intake Holes (per panel) 50 80 63 162 Height of Holes [m] 0.03 0.03 0.01 0.01 Width of Holes [m] 0.07 0.07 0.07 0.07 Hole area [m 2 ] 0.105 0.168 0.044 0.113 Hole area reduction due to jalousies (%) 55 55 0 0 Effective hole area [m 2 ] 0.047 0.076 0.044 0.113 Number of panels 1 3 4 4 Effective total volume of holes [m 2 ] 0.047 0.227 0.176 0.454 414 S. STANIŠIĆ, M. JEVTIĆ, B. DAS, Z. RADAKOVIĆ Table 2 Size and position of the modeled openings Inlet Cutout (C) Louvre side panel (B, A) Top Cutout (D) Area (width*height) [m 2 ] 0.282 x 0.632 0.36 x 0.602 0.101 x 0.632 x coordinate of center [m] 0.788 y coordinate of center (front side) [m] 0.012 y coordinate of center (rear side) [m] 1.307 z coordinate of center [m] 0.288 0.956 1.294 3. ASSUMPTIONS IN THE MODEL We experienced problems with the convergence and had to make simplifications to achieve model convergence. A final score is that 2 of 20 simulations converged with results in expected range, while 5 of 20 simulations were stopped due to very long computational time. The 2 successful simulations took 93 hours and 15 days, respectively. That is why we did not come to the point to include all relevant physical issues in the model. More precisely, FEM CFD model is focused on analyzing air flow and heat transfer due to the air mass transfer and heat transfer through the walls. Only the transformer compartment is modeled while it is assumed that the temperature in HV and LV compartments are equal to the ambient temperature. No thermal resistances to the heat conduction through the walls, through the ceiling and through the floor were considered. Kiosk floor is modeled as adiabatic. These approximations have smaller quantitative effect than two following approximations. The radiation heat transfer is not considered. Also, there is an approximation in the model for heat transfer along the radiators, which is important for the air buoyancy in the zone of the radiators. The basics about the buoyancy can be found in our previous publications [5] and [14]. Several improvements of the lumped model from [5] have been made in the meantime and will be published as an upgraded model for calculation of air buoyancy, based on the radiator modeled as a heat exchanger. In FEM CFD, fins are initially modeled as the rectangular aluminum bars. Transformer tank is modeled as a homogeneous body with low thermal conductivity (0.11 W / (mK)). This approximation causes the discrepancy of fin surface temperature from the real one, causing calculation error for the heat transfer to the air and calculation error for the air buoyancy. Two attempts were made at modeling air flow. The first was to consider laminar flow and the second with algebraic yPlus turbulent flow regime. No convergence with stationary solver has been achieved with either of the flow regimes, so transient solving was applied. The algebraic yPlus turbulence has been selected for the simulation of the flows inside closed areas [15]. It solves the flow everywhere and it is the most robust and least computationally intensive with good approximations for internal flow. Multiphysics comprising of heat transfer and turbulent flow, algebraic yPlus, was employed to the model. Effective size of kiosk openings has been modeled through boundary conditions in turbulent fluid flow physics: grille. This boundary condition incorporates effect of having a square mesh on the openings or louvres via head loss coefficient. Values for this coefficient are calculated using equations from [16]. FEM CFD Versus Lumped Thermal Model of Kiosk Substation with the Oil Immersed Distribution Transformer 415 After initial modeling of the fins as bars, they were reduced to surfaces, with thin layer boundary condition for heat transfer and interior walls for laminar flow. Comsol built-in default meshing sequence was used with the coarser setting. Stationary solver exhibited problems for both laminar and turbulent algebraic yPlus model. Online research and Comsol blog exploration pointed out to the experience that convective cooling sometimes poses a small transient that the stationary solver is not capable to solve. Transient (time domain) solver showed much better performance. After gathering such experience, we came to the idea to set initial value of the entire transformer block temperature to 74°C (equal to the steady-state top oil temperature, as presented in [5]), and thus to shorten the computational time needed to reach the steady-state (in respect to needed time if the initial condition would be the cold state). The initial temperature of the fins is set to 20°C. Algebraic yPlus turbulent model specific solver was used: transient with initialization. This solver is comprised of two study steps. The first step is to initialize the values of wall distances. This study step utilizes a fully coupled physics stationary solver in which initial values of all dependent variables (temperature, pressure, velocity field, wall distance in viscous units and reciprocal wall distance) are solved for using iterative GMRES method (Generalized Minimum RESidual). The second step uses segregated time dependent solver where first segregated node calculates for velocity field, pressure and temperature using iterative GMRES and the second node for wall distance in viscous units using direct PARDISO solver (PARallel DIrect Sparse Solver). Time range solved for is (0, 1, 60) [min], meaning that simulation is solved in minutes, starting from zero, with step 1 until one hour has been reached. Our assumption was that this period would be sufficient to reach a quasi-stationary state concerning air flow and cooling process since the entire domain volume of transformer already had the temperature as its steady-state condition. Nevertheless, because of setting the initial temperature of the fins to 20°C and adopting low thermal conductivity for the solid material of the tank no steady-state has been reached. 4. COMPUTING RESOURCES The best available computation resource we had was a single desktop with 64-bit Windows 7 Enterprise OS, Intel Core i5-6400 CPU and 32 Gb (2x16 Gb) of DDR4-2400/ PC4-19200 RAM. 5. CALCULATION RESULTS The following 6 Tables present the results of the post-processing of FEM CFD simulation results. Tables 3 and 4 show the differential pressures (differences of pressures of air exiting fins from above and air entering fins from bellow) averaged on the surface between each two fins, i.e. over/under the openings. The pressures on the transformer compartment openings in Table 4 are given in Pa as gauge pressures, i.e. the difference of absolute pressure and referent atmospheric pressure; referent pressure is 1.0133e5 Pa, at the level of the kiosk ceiling. 416 S. STANIŠIĆ, M. JEVTIĆ, B. DAS, Z. RADAKOVIĆ Table 5 shows the averaged temperatures on 10 surfaces leaning vertically on the fins and 2 surfaces leaning horizontally on the fins from below (z=0.508 m) and from above (z=1.308 m). Example of the bottom-most vertical surface on the front radiator is marked on Figure 2. Table 6 shows averaged temperatures on the openings. The ambient air temperature (outside the kiosk) is 20°C. Tables 7 and 8 present the results of the post-processing of FEM CFD simulation results for the total air flows in g/s through the surfaces as for averaged temperatures in table 5 (horizontal and vertical surfaces), i.e. cross-sections of the openings. Flow values are negative on surfaces where air is dominantly entering the radiator and positive on the outflow surfaces. Example of the bottom-most vertical surface on the front radiator is again the same as marked red on Figure 2. Table 3 Pressures difference (Pa) on the surfaces between bottoms and tops ΔP=Pexit-Pentry Between fins 1-2 2-3 3-4 4-5 5-6 6-7 7-8 Front side (8 fins) -9.3626 -9.3623 -9.361 -9.3599 -9.3601 -9.3612 -9.3646 Rear side (22 fins) -9.3603 -9.3619 -9.3623 -9.3624 -9.3617 -9.361 -9.3601 Between fins 8-9 9-10 10-11 11-12 12-13 13-14 14-15 Rear side (22 fins) -9.3584 -9.3543 -9.3537 -9.3575 -9.3617 -9.3637 -9.3642 Between fins 15-16 16-17 17-18 18-19 19-20 20-21 21-22 Rear side (22 fins) -9.365 -9.3656 -9.366 -9.366 -9.3644 -9.3622 -9.3579 Table 4 Pressures on the openings (in Pa as gauge pressures) Pressure [Pa] C B D Front (8 fins) 12.283 4.4086 0.4299 Rear (22 fins) 12.282 4.4089 0.4294 Table 5 An example of air temperature values [°C] over the fin height z-coordinate [m] Top surface 1.23 - 1.31 1.15 - 1.23 1.07 - 1.15 0.99 - 1.07 0.91 - 0.99 Front side (8 fins) 27.4 26.87 26.8 26.75 26.52 26.04 Rear side (22 fins) 26.82 26.11 25.87 25.77 25.67 25.39 z-coordinate [m] 0.83 - 0.91 0.75 - 0.83 0.67 - 0.75 0.59 - 0.67 0.51 - 0.59 Bottom surface Front side (8 fins) 25.63 25.19 24.75 24.51 24 23.66 Rear side (22 fins) 24.87 24.7 24.66 24.26 23.74 23.11 Table 6 The temperatures [°C] averaged on the openings C B D Front (8 fins) 19.84 21.77 28.67 Rear (22 fins) 19.91 21.69 28.39 FEM CFD Versus Lumped Thermal Model of Kiosk Substation with the Oil Immersed Distribution Transformer 417 Table 7 An example of air flow values [g/s] over the fin height z-coordinate [m] Top surface 1.23 - 1.31 1.15 - 1.23 1.07 - 1.15 0.99 - 1.07 0.91 - 0.99 Front side (8 fins) 2.6394 1.5214 1.2132 0.9277 0.5038 0.089 Rear side (22 fins) 5.8881 5.4824 4.2955 3.9995 3.1065 1.2586 z-coordinate [m] 0.83 - 0.91 0.75 - 0.83 0.67 - 0.75 0.59 - 0.67 0.51 - 0.59 Bottom surface Front side (8 fins) -0.2381 -0.4147 -0.4344 -0.4892 -0.428 -4.585 Rear side (22 fins) -0.3735 -1.1181 -0.6613 -0.1161 0.3325 -20.254 Table 8 The flows [g/s] on the openings Oppening C B D Front (8 fins) 42.84 14.413 40.188 Rear (22 fins) 69.67 18.666 39.236 Fig. 2 Example of the vertical surface (first out of 10 for front radiator) 418 S. STANIŠIĆ, M. JEVTIĆ, B. DAS, Z. RADAKOVIĆ Table 9 presents the values of the characteristic air flows. Table 9 Air flows Flows QC QHP QHTP QHNP QP QNP Values [g/s] 112.51 91.1 24.84 66.26 29.11 83.4 QC – Inlet opening C QHP – Upward flow around transformer through horizontal surface with coordinate z=0.51 m (just below the bottom edges of fins, QHTP – Upward flow through the horizontal surfaces below the fins (only the flow component entering both radiators from bellow) QHNP – Upward flow that does not enter the radiators: QHNP= QHP – QHTP QP – Total air flow into the radiators (QHTP is increased by the air entering from the side (see Table 7)) QNP – Part of flow through inlet openings (QC) minus flow entering into the radiators (QP): QNP= QC – QP Note: we suppose that the reason for the deviation of QHP from QC is the different mesh, i.e. the error caused by the interpolation for different meshes on the opening and on the horizontal plane below the fins. Large ratios QHP / QHTP and QNP / QP are the consequence of the air heating up on the tank surfaces which are not covered by the fins where air buoyancy also exists and friction is small (in fact, it appears only in velocity boundary layer. Total cooling surface in the lumped model is considered when the heat transfer coefficient (kp) is calculated. It is approximately supposed that the entire air mass, which is used for the calculation of the buoyancy and for the calculation of the frictional pressure drop in space between the radiator plates, flows vertically exclusively between the fins. Figure 3 presents the distribution of air velocity on the kiosk walls with the openings, being used to get the values in Table 8. Figure 4 visualizes space distribution of air velocity and temperature on the side with 8 fins; it is of relevance for values in Table 9. FEM CFD Versus Lumped Thermal Model of Kiosk Substation with the Oil Immersed Distribution Transformer 419 a) On the kiosk wall with the openings, side with 8 fins on the tank b) On the kiosk wall with the openings, side with 22 fins on the tank Fig. 3 Distribution of the air velocity (m/s) 420 S. STANIŠIĆ, M. JEVTIĆ, B. DAS, Z. RADAKOVIĆ Fig. 4 Space distribution of air flow pattern and temperature on the side with 8 fins FEM CFD Versus Lumped Thermal Model of Kiosk Substation with the Oil Immersed Distribution Transformer 421 6. THE TYPE OF THE RESULTS OBTAINED FROM FEM CFD AND LUMPED MODEL The characteristic temperatures, flows and the pressures can be obtained from the lumped model. The difference in respect to FEM CFD calculation is that lumped model delivers only one value for the pressure and the temperature at the bottom of the fins and also only one value at the top of the fins. The values from FEM CFD, which are to be compared with the ones from lumped model, are averaged on the surfaces. Since ideal upward air flow is supposed in the lumped model, there are no output values for air flow in and out through vertical surfaces in the zone of the fins (Tables 5 and 7). The lumped model results are presented in Tables 10 (temperatures), 11 (flows) and 12 (Pressures). Table 10 The values of temperatures obtained by lumped model Temperatures       Values [°C] 24.6 20.03 51.1 50.7 50.5 50 C  temperature on top of inner side of opening C BR  temperature on entry to the radiator TR  temperature on exit from the radiator ceil  temperature on kiosk ceiling near the kiosk wall D  temperature on top of opening D B  temperature on top of opening B Table 11 The values of flows obtained by lumped model Flows QC QBR QD QB QKDB QKBC Values [m 3 /h] 0.0322 0.129 0.0226 0.01 0.0205 0.0004 There are 4 openings C, 4 openings D and 4 openings B (flow through opening A is practically the same as through B, so it is considered as there is 4 openings B instead of 3 B and 1 A) QC  flow through opening C QBR  flow through the radiator QD  flow through opening D QB  flow through opening B QKDB  flow downstream the kiosk wall (between the openings D and B) QKBC  flow downstream the kiosk wall (between the openings B and C) Table 12 The values of pressure differences obtained by lumped model Press. diff.      Values [Pa] -2.662 -8.5479 -1.7405 1.7434 0.4671 Press. diff.     Values [Pa] 0.2726 2.2545* (2.486)** 0.0675 9.1455* (9.3594)** * on the inner side of the kiosk, ** on the outer side of the kiosk 422 S. STANIŠIĆ, M. JEVTIĆ, B. DAS, Z. RADAKOVIĆ pC-BR  pressure difference between middle of opening C  entry to the radiator pBR-TR  pressure difference between entry to the radiator  exit from the radiator pTR-Ceil  pressure difference between exit from the radiator  kiosk ceiling pCeil-D  pressure difference between kiosk ceiling  middle of opening D pDIn-DOut  pressure difference between inner side of middle of opening D  outer side of middle of opening D pBIn-BOut  pressure difference between inner side of middle of opening B  outer side of middle of opening B pD-B  pressure difference between middle of opening D  middle of opening B pCOut-CIn  pressure difference between outer side of middle of opening C - inner side of middle of opening C pB-C  pressure difference between middle of opening D  middle of opening C 7. CONCLUSIONS FEM CFD method is relatively new and presents a powerful tool for analysis of wide variety of heat transfer problems including fluid flow. Nevertheless, as presented in the paper, severe convergence problems can appear when using software based on this method. From that point of view, publishing the practical experience of its application is valuable. Convergence problems we encountered were clearly stated in the paper. At the end, FEM CFD results that were obtained only gave us qualitative representation of air flow. It was not possible to compare the results with the results of lumped model since the output quantities were not the same. In the experiment there was no record which corresponds to the FEM CFD simulation (initial state is different and no steady-state has been reached in FEM CFD). Stronger computational resources could probably make it feasible to use smaller mesh and to increase convergence. Another option, combined with the previous one, is to perform custom meshing in critical zones, i.e. not to use automatic mesh generation as we did. Such work is presented in [17], where similar problem, but for transformer placed under the ground surface, is considered using FEM CFD. The approach in [17] is stricter with fewer simplifications, but with much stronger hardware recourses solving a model with much higher mesh cells number (as well as performing grid independence verification). At the end, the following conclusions about air flow distribution were drawn from the simulations that converged, which could not have been seen in lumped model developed and applied in our previous work: 1. A part of the air exits the radiator before it reaches the top of the radiator, streaming toward the outlet cutouts. Similar situation happens for the air entry: part of air flows from the inlet openings and enters the radiator on the vertical boundary surface of the air ducts between the fins. 2. There is significant upward air flow outside the zone of the fins, caused by the buoyancy in the areas of the tank surfaces which are not covered by the fins (see Section 5). These finds should be kept in mind while building the lumped models, i.e. a way of considering their influence (via elements of lumped model) should be explored. FEM CFD Versus Lumped Thermal Model of Kiosk Substation with the Oil Immersed Distribution Transformer 423 REFERENCES [1] K. Baral and I. Primus, "Lebensdauer eines 630 kVATransformators in einer Beton-Netzstation," Elektrizitaetswirtschaft, vol. 8, pp. 268–276, 1979. [2] I. Primus, "Temperaturen in Netzstationen–wirtschafliche Bedeutung und Einfluss Factoren," Elektrizitaetswirtschaft, vol. 16, pp. 451–460, 1976. [3] I. Primus, "Temperaturen in Netzstationen–Messergebnisse und Deutung," Еlektrizitaetswirtschaft, vol. 22, pp. 833–842, 1976. [4] Z. Radakovic and S. Maksimovic, "Non–stationary thermal model of indoor transformer stations," Electrical Engineering (Archiv fur Elektrotechnik), vol. 84, no. 2, pp. 109-117, 2002. [5] Z. Radakovic, M. Jevtic and B. Das, "Dynamic thermal model of kiosk oil immersed transformers based on the thermal buoyancy driven air flow," International Journal of Electrical Power & Energy Systems, vol. 92, pp. 14-24, Nov. 2017. [6] A. K. Das and S. Chatterjee, "Finite element method-based modelling of flow rate and temperature distribution in an oil-filled disc-type winding transformer using COMSOL multiphysics," IET Electric Power Applications, vol. 11, no. 4, pp. 664 - 673, 2017. [7] Y. Huijuan, Y. Tingfang, X. Rui and P. Chunhua, "Numerical Simulation of Ventilation for Main Transformer Room of Indoor Substations," The Open Automation and Control Systems Journal, vol. 7, pp. 630-639, 2015. [8] H. Liu, Y. Hao, M. Fu, D. Wang and L. Yang, "Study on ventilation of indoor substation main transformer room based on COMSOL software," in Proceedings of the 1st International Conference on Electrical Materials and Power Equipment (ICEMPE), Xi'an, China, July 2017. [9] M. R. Nalamwar, D. K. Parbat and D. Singh, "Study of effect of windows location on ventilation by CFD simulation," International Journal of Civil Engineering and Technology, vol. 8, pp. 521-531, 2017. [10] T. Yu, H. Yang, R. Xu and C. Peng, "Simulation Study on Ventilation & Cooling for Main Transformer Room of an Indoor Substation," Journal of Multimedia, vol. 9, no. 8, pp. 1040-1047, 2014. [11] M. Banjac, "Application of Computational Fluid Dynamics in cooling systems design for special purpose objects," FME Transaction, vol. 42, pp. 26-33, 2014. [12] S. B. Paramane, W. V. d. Veken and A. Sharmac, "A coupled internal–external flow and conjugate heat transfer simulations and experiments on radiators of a transformer," Applied Thermal Engineering, vol. 103, pp. 961-970, June 2016. [13] S. B. Paramane, K. Joshi, W. V. d. Veken and A. Sharma, "CFD Study on Thermal Performance of Radiators in a Power Transformer: Effect of Blowing Direction and Offset of Fans," IEEE Transactions on Power Delivery, vol. 29, no. 6, pp. 2596-2604, Dec. 2014. [14] Z. Radakovic and M. Sorgic, "Basics of Detailed Thermal-Hydraulic Model for Thermal Design of Oil Power Transformers," IEEE Trans. on Power Delivery, vol. 25, no. 2, pp. 790-802, 2010. [15] https://www.comsol.com/blogs/which-turbulence-model-should-choose-cfd-application/. [16] I.E. Idelchick, Handbook of hydraulic resistance, 3rd ed., Florida: CRC Press Inc., 1994. [17] J.C. Ramos, M. Beiza, J. Gastelurrutia, A. Rivas, R. Anton, G. S. Larraona and I. de Miguel, "Numerical modelling of the natural ventilation of underground transformer substations," Applied Thermal Engineering, vol. 51, no. 1-2, pp. 852-863, March 2013.