Thermal management of a distribution transformer: An optimization study of the cooling system using CFD and response surface methodology

Thermal management of a distribution transformer: An optimization study of the cooling system using CFD and response surface methodology

Electrical Power and Energy Systems 104 (2019) 443–455 Contents lists available at ScienceDirect Electrical Power and Energy Systems journal homepag...

3MB Sizes 0 Downloads 37 Views

Electrical Power and Energy Systems 104 (2019) 443–455

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Thermal management of a distribution transformer: An optimization study of the cooling system using CFD and response surface methodology

T



Leyla Raeisiana, Hamid Niazmanda, , Ehsan Ebrahimnia-Bajestanb, Peter Werlec a

Department of Mechanical Engineering, Faculty of Engineering, Ferdowsi University of Mashhad, Mashhad, Iran Department of Mechanical Engineering, Quchan University of Technology, Quchan, Iran c Department of Electrical Power Systems, Schering-Institute, Leibniz University of Hannover, Hannover, Germany b

A R T I C LE I N FO

A B S T R A C T

Keywords: Transformer Thermal management Hotspot Optimization Response surface method

In this paper, a numerical scheme has been developed to examine the effective parameters on thermal management of distribution transformers and subsequently to optimize their cooling systems. In this regard, the response surface methodology (RSM) was used as the optimization method to minimize the hotspot temperature in the transformer as the response factor. A comprehensive three-dimensional computational scheme was employed considering the detailed geometrical specifications of an actual 200kVA distribution transformer to obtain the temperature field and the hotspot temperature. The accuracy of the numerical model was established via comparing the numerical results with the measured temperatures of a running transformer. The thermal variations of the thermo-physical properties of the transformer oil are determined experimentally and incorporated in the numerical modeling. A comprehensive parametric study among seven evidently effective parameters has been performed to identify the most effective parameters on the thermal performance of the transformer. It was found that fin height, length, and spacing are the more influential parameters among the examined parameters, which are also considered as the input variables in the optimization procedure. According to the RSM, the effects of the variations of these input variables in pre-specified ranges on the response, which is the hotspot temperature, are examined through the suggested runs by RSM. The results indicated that the hotspot temperature is more influenced by the fin height as compared to the fin length and spacing. Furthermore, the hotspot temperature decreases with the increase in fin height and length, while decreases as fin spacing increases. In addition, a correlation for the variations of the hotspot temperature as a function of the fin height (H), length (L), and spacing (S) is suggested using RSM. The significant finding is that the proposed optimum transformer configuration (H = 0.9 m, L = 0.08 m and S = 0.036 m) leads to the hotspot temperature reduction of about 16 °C as compared to the actual transformer geometry, which greatly affects the transformer life expectancy and safer performance.

1. Introduction Electric transformers as an essential electric equipment perform a vital function in supplying the electric energy at required voltage levels to consumers. An under-loaded transformer not only meets an electrical process but also experiences a thermal activity that is impelled by heat generated due to the losses in core and windings. The generated heat is the major cause of transformer temperature rise and its aging [1] and should be removed by employing a suitable cooling system. Among the various parameters that affect the transformer performance, the location of the maximum temperature of the solid insulation named ‘hotspot temperature’ (HST) and its temperature level have been identified as the main reasons of the transformer aging and failures [2]. In the



hotspot location the depolymerization process of cellulose insulation gradually occurs, which degrades the mechanical properties of cellulose paper such as tensile strength and elasticity, which in turn, makes the paper brittle and incapable of enduring electrical forces and vibrations. This irreversible deterioration strongly affects the transformer end of life [3]. Consequently, thermal management of transformers by means of appropriate cooling systems is the key factor of their safe operations and improving their productivities [4]. The thermal management efficiency of a transformer relies on various parameters such as dielectric strength and thermal properties of insulating media, operational conditions and geometrical design [5]. Among the mentioned factors that influence the hotspot and life expectancy of transformer, system configuration is the important one. The design of core and windings is

Corresponding author. E-mail address: [email protected] (H. Niazmand).

https://doi.org/10.1016/j.ijepes.2018.07.043 Received 14 January 2018; Received in revised form 30 May 2018; Accepted 19 July 2018 0142-0615/ © 2018 Elsevier Ltd. All rights reserved.

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Nomenclature As Cp Ch CH CL CW FT g h H k LC L LVW HVW Nu p Ps Ra S T TH

TL TW u v w WH W Y1 Y2

cross sectional area (m2) specific heat capacity (J.kg-1.K−1) oil channel core height core length core square side fin thickness gravitational acceleration (m.s−2) convective heat transfer coefficient (W.m−2.K−1) fin height (m) thermal conductivity (W.m−1.K−1) characteristic length (m) fin length (m) low voltage windings high voltage windings Nusselt number pressure (Pa) perimeter (m) Rayleigh number fin spacing(m) temperature (K) tank height (m)

tank length(m) tank width (m) velocity in X direction (m.s−1) velocity in Y direction (m.s−1) velocity in Z direction (m.s−1) winding height (m) oil channel width (m) active part vertical position (m) fin vertical position (m)

Subscripts a s ∞

air layer surface ambient

Greek letters β ρ μ α ν

thermal expansion coefficient density (kg m−3) viscosity (Pa s) thermal diffusivity (m2 s−1) Kinematic viscosity (m2 s−1)

and Salama [12], presented a method of estimating the lifetime of transformer insulation according to some artificial histories for the ambient temperature and the transformer load, which were used as the inputs to the Monte Carlo simulation technique in order to find the thermal lifetime of the insulation of a given transformer. Although, estimation of maximum temperature and lifetime of transformers are of great practical importance, simplifications of the governing equations and geometrical models in above-mentioned studies lead to the case sensitive and relatively low accurate results. Literature review shows that most of the previous studies are belong to the second group, where the impacts of geometry and coolant flow inside the channels of the windings are considered on the temperature distribution of windings, using a 2D or an axisymmetric model. Clearly, in these models, which are mostly developed for power transformers, the geometry of windings are basically examined not the entire transformer. In the following, some of the papers related to this category are briefly reviewed. Hosseini et al. [13] investigated the effects of oil channels thickness and washer arrangements on the hotspot and also temperature distribution in windings using the FLUENT software. Besides, Torriano et al. [14] numerically studied the impacts of the operational parameters, such as the mass flow rate and boundary conditions on the flow and temperature distributions of the windings. They concluded that in higher flow rate, the better cooling of the windings is obtained, since the windings temperature distributions become more uniform and the hotspot temperature reduces. Similarly, Zhang and Li [15] demonstrated that the coolant flow path arrangement plays a vital role on the cooling of winding disks, such that the maximum temperature usually occurs in the disks surrounded by the ducts with minimum coolant flow rate. Toriano et al. [16] also performed a three dimensional (3D) conjugate heat transfer simulation of oil flow through the cooling ducts and demonstrated the importance of 3D modeling. They reported the weakness of the axisymmetric modeling due to the missing of annular velocity and temperature gradients in oil channels of windings. Campelo et al. [17] applied Computational Fluid Dynamic (CFD) and Thermal-Hydraulic Network Models (THNMs) to study and compare the flow and temperature distributions in the 2D and 3D modes. Skillen et al. [18] investigated the oil flow and temperature distributions in the low voltage winding of a power transformer using the CFD approach along with the temperature dependent coolant

usually performed by two main structures, which are the layers model for distribution transformers and disks for power ones, while reviewing the literature indicates that various designs can be used for cooling systems. A number of researchers have investigated the effects of cooling system design on the overall performance of the transformers. These efforts can be categorized according to the methodology of the investigation into three main groups as: (i) estimation of the hotspot temperature and location as well as the lifetime of transformers; (ii) study of temperature distribution in windings, and (iii) simulation of convective heat transfer of insulating media in transformers. In the first group of studies, simplified forms of energy equation are usually solved along with the proper boundary and operational conditions, where the complications of the realistic geometry are ignored. This method of study provides the possibility of fast prediction of the hotspot. For instance, Taghikhani and Gholami [6] investigated the location and value of a power transformer's hotspot by solving the partial differential form of heat conduction equations numerically for a winding disk transformer, which were well compatible with the results of the actual tests on site. Pradhan and Ramo [7] also presented a theoretical model to examine the hotspot, based on a boundary value problem of heat conduction in windings employing the finite integral transform technique. Susa and Nordman [8] presented a simple model, based on the exponential iterative calculation procedure, for calculating the hotspot temperature by relating the hotspot temperature to the top oil temperature. They also compared their results with the IEEE Annex G method, where reasonable agreements were observed. Furthermore, Hajidavalloo et al. [9] assessed the effect of sun shield on the hotspot temperature of a transformer, both experimentally and numerically. It was observed that the oil temperature under the shadow mode was remarkably lower than that of the normal mode, as the reduction of hotspot temperature in summer was as high as 7 °C and the lifetime enhancement by the sun shield was nearly 24%. Souza et al. [10] applied an evolving fuzzy model to study the lifetime of distribution transformers. Koufakis et al. [11] developed a theoretical method to calculate the insulating resistance of the insulation paper in distribution transformers at several temperatures, which strongly influenced the transformer lifetime. They also presented a method to calculate the transformer life cycle. In a similar research, Abu-Alanian 444

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Gastelurrutia et al. [29], they investigated the thermal behavior of a distribution transformer by a turbulence model considering the oil flow in a thin slice of transformer. In the third category, more detailed and accurate results are expected to be obtained, however, some simplifications in the geometry or the operational conditions may limit the accuracy of the models. Surveying the literature indicates that, although a reasonable amount of work has been carried out on transformer's cooling, investigating the effects of various geometrical parameters on the hotspot temperature and optimizing the cooling system are still in its initial stages. To the best knowledge of the authors, a 3D simulation of convective heat transfer that considers the complete geometry of the transformer such as hollow fins and details of oil channels to examine the effects of the geometrical parameters on the oil thermal behavior and the optimization of the transformer geometry has not been performed previously. Accordingly, the main objective of the present study is to introduce an optimum geometry from the thermal point of view aiming at lowering the hotspot temperature in comparison with the conventional design. This target was achieved by the following steps. Firstly a 3-D numerical simulation of the oil thermal behavior inside a 200 kVA distribution transformer has been performed by the ANSYS Fluent software considering all details in the fluid domain. The numerical model was also validated with the experimental data. Then, the effects of different geometrical parameters were numerically investigated in order to identify the most influential parameters on the hotspot temperature. Finally, an optimized geometry was generated for the transformer based on the recognized effective parameters using the response surface method and the design expert software.

properties. They mainly showed that the CFD approach is capable of capturing the oil flow and temperature patterns, more accurately than the common industrial methods, in which the flow velocity is assumed to be uniform. Codde et al. [19] indicated that the common THNMs are not capable of accurately predicting the hydrodynamic flow properties in zig-zag cooled transformer windings. They extracted several correlations from a CFD simulation to improve the accuracy of the network models for the zig-zag type of transformer's windings. Yatseveski et al. [20] also analyzed the temperature distributions of the windings by FLUNET software and compared them with a running transformer. In another work, Rahimpour et al. [21] also studied the temperature distribution dependency on the washer design in a disk type transformer. Furthermore, Smolka [22] by employing the CFD and genetic algorithm presented an optimum configuration for the windings and their oil channels in order to decrease the hotspot temperature. Consequently, the study of the second category demonstrated that simulation of windings as solid disk provides significant information for designing the efficient configurations for the oil channels located inside the windings, yet, the analysis of the complete cooling system and investigating the importance of effective factors in a full 3D transformer considering core and windings is inaccessible. In addition, employing simple thermal and hydrodynamic boundary conditions for simulation of the windings are the main source of errors in the modeling of the hotspot temperature and location in transformers. In order to overcome the abovementioned limitations, the researchers have been encouraged to model the complete convective heat transfer of the entire transformer in the ambient conditions. However, due to the complexity of this method, which falls into the third group, limited studies are available as reviewed in the followings. El Wakil et al. [23] employed a numerical method to study the heat transfer and fluid flow of an insulating media through the one phase of a 3-phase 40 MVA power transformer. They considered six different geometrical configurations along with the six different oil flow rates in order to obtain an optimum cooling system. Oh and Ha [24] considered an axisymmetric numerical model for turbulent natural convection in a 20 kVA one phase transformer. Radakovic et al. [25] investigated some parameters like wind velocity on the oil temperature near the outer surface of a kiosk transformer. Furthermore, Tsilia [26] conducted a numerical study using a 3D finite element method to predict the temperature field of the transformer. Chereches et al. [27] performed a numerical parametric study on a step down power transformer to investigate the effects of various parameters on the performance of the cooling system. They reported that placing an obstacle at the bottom of the transformer helps the coolant to flow towards the hottest part and decreases the maximum temperature. Eslamian et al. [28] simulated the heat transfer of dry typed transformer by FEM as well as CFD models and investigated the effects of several factors such as duct width on temperature distribution of transformers. In another research by

2. Numerical modeling 2.1. Device In this work, a well-known 200 kVA distribution transformer was considered as the case study, which was manufactured in 2015 by Iran Transfo Company. The studied transformer has an active part consist of the core and windings, as well as a cooling part including oil, oil channels, tank and hollow fins attached to the tank as can be seen in Figs. 1a–4a. Furthermore, its low voltage (LV) and high voltage (HV) windings are 400 V and 20 kV, respectively. 2.2. Domain In the considered transformer, mineral oil circulates through the active part and tank as a cooling media in order to absorb heat from the core and windings and rejects it to the surrounding air via tank and fin surfaces, to provide a reliable level of temperature distribution in the

Fig. 1. Specification of the transformer tank, a) the real transformer, and b) the simulated geometry of the tank and hollow fins. 445

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

calculation, which is relatively far from 109 related to the start of turbulent flows. In the rest of the transformer cavity, the flow velocity is rather low and the oil passages are narrow enough to satisfy the laminar flow condition. It is also assumed that the oil flow is steady state, incompressible and three dimensional with temperature dependent thermo-physical properties, based on the operational conditions of oil inside the transformer. The generated heat in the core and windings are modeled as a constant heat flux at their interfaces with the transformer oil as will be discussed in more details in the boundary conditions section. Considering all these assumptions the governing equations for the oil are as follows. Continuity:

transformer. In this regard, the oil inside the transformer was considered as the numerical domain and its convective heat transfer was three dimensionally modeled and simulated. 2.3. Geometry In the present study, thermal behavior of the oil flow inside the transformer has been simulated according to the geometrical and operational conditions of the mentioned real transformer. The body of the actual transformer and the corresponding simulated geometry have been shown in Fig. 1. Furthermore, the internal details of the transformer including all oil channels inside the active part and hollow fins with their exact dimensions were considered, according to the Table 1, Figs. 2 and 3. Considering Fig. 2a, the cross section of the core in the actual transformer is more like a polygon, however, for simplification, it was modeled as a square [26]. In Fig. 2, all oil channels are shown, including the narrow oil channel located between the core and LV winding (Ch. 1), the one in the LV winding (Ch. 2), and the channel between LV and HV windings (Ch. 3), as well as the HV winding oil channel (Ch. 4). These oil passages are fabricated for appropriate heat removal from the core and windings. Additionally, Fig. 4 displays the hollow fins and their corresponding modeled geometry.

∂ (ρw ) ∂ (ρv ) ∂ (ρu) =0 + + ∂z ∂y ∂x

(1)

X direction Momentum:

∂ (ρu2) ∂ (ρuv ) ∂ (ρuw ) ∂P ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ + + =− μ + + μ ∂x ∂y ∂z ∂x ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂ ⎛ ∂u ⎞ + μ ∂z ⎝ ∂z ⎠ ⎜

2.4. Mesh



(2)

Y direction Momentum:

Mesh generation of the computational domain was carried out in ANSYS ICEM, which provides different types of meshes. For the considered geometry, however, both structured mesh and unstructured mesh are possible. Yet, according to the authorś numerical experience on the simpler transformer geometries, the unstructured meshes are not only easier and faster in generation in comparison with the structured mesh, but also they provide the same accuracy as the structured mesh, and therefore, it was employed for meshing the fluid domain. In this regard, according to the geometry, the generation of mesh was implemented in two steps. Firstly, the tetragonal (TET) elements are generated throughout the whole geometry and then, by using the “TET to HEX conversion” option, the possible TET elements were converted to the hexagonal (HEX) elements in order to decrease the mesh numbers, and therefore, the running time. Consequently, in relatively larger areas the TET elements were converted to HEX elements, while in narrow regions like fins and oil channels TET elements were retained as shown in Fig. 5. Then, in order to find the appropriate mesh size, a grid independency study was carried out on the values of the average and maximum temperature of the transformer, as presented in Table 2. In this grid study, the number of cells for every mesh system is about 1.5 times more than that of the previous mesh system. It can be seen that the change in the maximum temperature from Grid 1 to Grid 2 system is about 5 °C, however, the maximum temperature remains nearly constant from Grid 2 to Grid 4 system. The average temperature also gradually reduces from Grid 1 to Grid 4, with almost the same values for Grid 3 and Grid 4. Considering the accuracy and the computational cost of the numerical solution, the mesh system of Grid 2 with the total cell number of 2,948,836 was selected as the mesh system for all numerical simulations of this study.

∂ (ρvw ) ∂ (ρv 2) ∂ (ρuv ) ∂P ∂ ⎛ ∂v ⎞ ∂ ⎛ ∂v ⎞ =− + + + + μ μ ∂z ∂y ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂y ∂x ∂ ⎛ ∂v ⎞ + μ −ρg ∂z ⎝ ∂z ⎠ ⎜



(3)

Z direction Momentum:

∂ (ρw 2) ∂ (ρvw ) ∂ (ρuw ) ∂P ∂ ⎛ ∂w ⎞ ∂ ⎛ ∂w ⎞ =− + + + + μ μ ∂z ∂z ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂y ∂x ∂ ⎛ ∂w ⎞ + μ ∂z ⎝ ∂z ⎠ ⎜



(4)

Energy: Table 1 Geometrical specifications of the base transformer.

2.5. Governing equations In the present problem, natural convective heat transfer is the predominant mechanism of transferring heat from electric parts to the oil, and also from outer surfaces of the transformer to the ambient air. Furthermore, according to the transformer geometry and the calculated Rayleigh number, the oil flow in the transformer was considered as laminar flow. It should be mentioned that the critical part of the transformer which has the potential to be in the turbulence regime is the top section that is the gap between the core and top wall, which has a maximum Rayleigh number of about Ra = 5 × 108 according to the 446

Parameters

Value (m)

Fin length (L) Fin thickness (FT) Fin spacing (S) Fin height (H) Tank length (TL) Tank width (TW) Tank height (TH) Core height (CH) Core width (square side) (CW) Core length (CL) Windings height (WH) Inner diameter of LVW1 Outer diameter of LVW1 Inner diameter of LVW2 Outer diameter of LVW2 Inner diameter of HVW1 Outer diameter of HVW1 Inner diameter of HVW2 Outer diameter of HVW2 Average distance between the core and LVW1; oil channel (Ch. 1) Oil channel width (Ch. 2) Oil channel width (Ch. 3) Oil channel width (Ch. 4) Vertical distance of the active part from the upper surface of tank (Y1) Vertical distance of the fin from the upper surface of tank (Y2)

0.08 0.007 0.04 0.7 1.02 0.42 0.90 0.68 0.105 0.68 0.36 0.078 0.0855 0.0895 0.097 0.107 0.119 0.124 0.138 0.015 0.005 0.01 0.007 0.18 0.03

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Fig. 2. Specifications of the active part and its oil channels, a) the real transformer b) the simulated geometry.

Fig. 3. 3D specifications of the core, windings and oil channels, a) the real transformer, b) the simulated geometry, and c) positioning of the active part and fins.

∂ (ρCp uT ) ∂x

+

∂ (ρCp vT ) ∂y

+

∂ (ρCp wT ) ∂z

∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ k + k ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂ ⎛ ∂T ⎞ + k ∂z ⎝ ∂z ⎠

=



to 100 °C by DSC1 device made by Mettlertolendo, under ASTME 1269 standard. To examine the thermal conductivity, a KD2 Pro made by Decagon Devices Inc. was employed to measure thermal conductivity at temperatures of 25, 35 and 45 °C, according to ASTM D5334-14. Besides, the viscosity of the oil was evaluated at 20, 40, 60 and 80 °C and constant shear rate of 100 (1/S) by Brookfield DV-III Ultra Programmable Rheometer. The obtained correlations based on the measured data, which have been implemented in the numerical modeling, are presented in Eqs. (6)–(9), where the temperature, T is based on the Celsius unit.



(5)

2.6. Properties of the coolant oil The flow and temperature distributions are strongly dependent on the oil thermo-physical properties. Therefore, it is important to exactly determine the thermal variations of the properties, which are commonly ignored in the previous studies. For this purpose, the density measurements were carried out using an Anton Paar DMA 5000 oscillatory densitometer at temperatures in the range of 30 °C–80 °C, according to ISO 17025 standard. Also, specific heat capacity at constant ambient pressure was measured for temperatures in the range of 20 °C 447

ρ (T ) = −0.6881T + 897.62

(6)

Cp (T ) = 6.6096T + 1600.6

(7)

K (T ) = 0.0009T + 0.1035

(8)

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

2.7. Boundary condition In transformers, the electrical and magnetic losses of the active part including the hysteresis and eddy current losses of the core and copper loss of the winding are the main sources of heat generation. The generated heat can be obtained either by solving the electromagnetic field of the core and windings, or through the short circuit and open circuit tests of the core and windings. The solution of the electromagnetic field provides the heat generation distribution, which appears as a source term in the solution of the transformer thermal field, which must be performed in an iterative manner. Considering the fact that the core and windings are made of iron and copper with high thermal conductivities, it is reasonable to assume that the generated heat losses are distributed uniformly within their volume. Therefore, the solution of the electromagnetic field is not required for the solution of the transformer thermal field, since the short circuit and open circuit tests can provide the power losses of the core and windings, respectively. However, the solution of the heat conduction inside the core and windings, similar to the previous case, is still required its interface temperature with the oil as the boundary condition. Clearly, the solution of this conjugate problem is also time consuming due to the unknown interface temperature, which must be determined through an iterative process. The computational cost associated with the above mentioned thermal modeling systems, when it is combined with an optimization process was practically beyond the capabilities of the computing facilities available to the present work. Considering the aim of the present study, which is the optimization of the transformer geometry by the response surface method, assuming the generated heat losses as a constant and uniform heat flux boundary condition at the interface of the oil with the active parts is found to be a remedy that provides the adequate level of accuracy with a reasonable computational cost. In this regard, under the full load operating condition of the studied transformer, the heat fluxes for the windings and core were considered as 470 and 380 W.m−2, respectively. Furthermore, based on the realistic ambient conditions of the transformer used for validating the numerical model, natural convection is the mechanism of heat transfer from the transformer surfaces to the ambient air. The details of the surface boundary conditions and related equations are taken from [30] and [31] and have been summarized in Table 3, according to the surface geometry and its orientation. In this table, Lc is the characteristic length of the examined surface, which is used in the Ra and Nu correlations. The thermo-physical properties of air at two ambient temperatures considered in the present study are listed in Table 4. In Table 3, the proposed equation for the Nusselt number of the finned surfaces is valid for the case of natural convection cooling of vertical finned surfaces. In this semi-empirical equation, the first term in right-hand side of the equation corresponds to the fully developed natural convection between the parallel plates with small gaps (and relatively long channels), however, the second term defines the natural convection heat transfer from an isolated vertical plate (valid for the inlet region and in relatively short channels), which becomes the dominant term for the relatively larger fin gaps. These two different contributions were combined by Bar-Cohen and Rohsenow [32] using the correlating expression suggested by Churchill and Usagi [33], which is valid for a wide range of vertical finned surfaces with various gaps.

Fig. 4. 3D specifications of the hollow fins, a) the real transformer, and b) the simulated geometry.

Fig. 5. Mesh distribution at the Z = 0 cross section plane.

Table 2 Grid independency tests for the full load transformer. Item

Value

Grid system Maximum mesh size(mm) Cell numbers Maximum temperature (K) Average temperature (K)

Grid 1 0.01 1,718,407 383.82 372.63

1055 μ (T ) = exp ⎛ −3.591⎞ ⎝ T + 141.8 ⎠

1000

Grid 2 0.008 2,948,836 378.78 372.14

Grid 3 0.007 4,108,029 379.312 371.57

Grid 4 0.006 6,196,590 379.568 371.49

2.8. Numerical implementation The control volume based commercial software of ANSYS FLUENT was employed for the solution of the mentioned governing equations, with the prescribed boundary conditions. Discretizations of the governing equations of momentum and energy were carried out by the second order upwind scheme. Presto and SIMPLE were used for the pressure and pressure–velocity coupling equations, respectively [34].

(9)

448

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Table 3 The transformer surface boundary conditions with the ambient air at temperature of T∞ =28 °C. Surface

Nusselt correlation

L

Upper Bottom Vertical finned

Nu = 0.59Ra0.25[30] Nu = 0.27Ra0.25[30]

As/Ps As/Ps S

Nu = [

576 (Ra LC / H )2

+

2.873 ]−0.5 (Ra LC / H )0.5

[30]

Table 4 Air properties at the standard pressure of 1 atm [30]

C,

definition

L C(m), value

Calculated h (W.m−2.K−1)

0.149 0.149 0.04

6.25 3.00 4.45

Table 5 Boundary conditions for the validation cases of the present model.

T (°C)

β (K−1)

ν (m2.s−1)

α (m2.s−1)

k (W.m−1.K−1)

15 28

0.0035 0.0033

1.47 × 10−5 1.585 × 10−5

2.009 × 10−5 2.174 × 10−5

0.02476 0.02570

Surface

Values related to the boundary conditions

Upper horizontal Bottom horizontal Vertical finned Windings Core

In addition, all the under relaxation factors are gradually increased from the initial value of 0.01 to the final value of 0.45 to avoid any divergence during the calculations, which is one of the major challenges in the simulation of such natural convection problems. In the postprocessing step, the results were analyzed to obtain the hotspot temperature of the transformer. It is also worth to be noted that the computational cost of each run was about 156 h, employing a desktop computer with specifications of Intel Core(TM) i7-4790 K CPU with 32 GB RAM.

10% full load T∞ = 15 °C

30% full load T∞ = 21 °C

h = 4.3 W m−2 K−1; h = 2.1 W m−2 K−1; h = 3.1 W m−2 K−1; q = 4.5 W m−2 q = 380 W m−2

h = 4.8 W m−2 K−1; h = 2.23 W m−2 K−1; h = 3.27 W m−2 K−1; q = 46.2 W m−2 q = 380 W m−2

Table 6a The temperature of specified points on the surface of the transformer predicted by the present numerical model as compared to the measured data at 10% full load.

2.9. Validation In order to validate the current numerical model, the transformer surface temperatures were compared with the measured data obtained from an actual 200 kVA distribution transformer operating at 10% and 30% of its full load. For this purpose, a thermography camera was used to record the temperature of the outer surface of the transformer. For the transformer at 10% of its full load, 12 points along a horizontal line on the surface of the transformer tank with the constant intervals of 92.7 mm as shown in Fig. 6a were selected as the measuring points. However, for the case of 30% of the full load, 10 points along a vertical line on the edged fin surface, as shown in Fig. 6b, were chosen with the constant intervals of 77.8 mm. The measurements were performed in Mashhad city of Iran at the end of September at 2 pm for 10% of full load and in Kerman in the middle of August at 6 am for 30% of full load, when the transformers were in shadow and the wind speed was almost negligible. Considering the operating and ambient conditions during the measurements, the employed boundary conditions in numerical solution are listed in

Location

Measured data (K)

Numerical results (K)

Deviation (K)

1 2 3 4 5 6 7 8 9 10 11 12

301.05 302.25 302.55 302.25 302.05 302.15 302.15 302.35 301.55 301.55 301.15 301.15

302.945 302.945 302.094 302.945 302.945 302.945 302.945 302.945 302.094 302.945 302.945 302.945

1.895 0.695 −0.456 0.695 0.895 0.795 0.795 0.595 0.544 1.395 1.795 1.795

Table 5. In Tables 6a and 6b the measured data and the calculated temperatures are listed for both cases. Reasonable agreements between the results of both cases indicate the accuracy of the present numerical model.

Fig. 6. Locations of the temperature measured points for validation of the numerical solution a) 10% full load b) 30% full load. 449

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Table 6b The temperature of specified points on the surface of the transformer predicted by the present numerical model as compared to the measured data at 30% full load. Location

Measured data (K)

Numerical results (K)

Deviation (K)

1 2 3 4 5 6 7 8 9 10

315.05 315.65 314.95 313.55 314.35 313.85 313.25 312.9 311.35 311.05

316.137 317.485 317.485 317.035 316.812 316.586 316.362 315.912 315.463 313.665

1.087 1.835 2.535 3.485 2.462 2.736 3.112 3.012 4.113 2.615

Table 7 Effects of geometrical variations on the hotspot temperature (HST). Parameters and values

HST (K)

HST variations with respect to the base case

Fin length, L (m) L = 0.08 L = 0.1 L = 0.12

386.5 376.8 370.2

Base case −9.6 −16.4

Fin thickness, FT (m) FT = 0.007 FT = 0.008 TF = 0.009

386.5 384.3 383.7

Base case −1.7 −2.4

386.5 386.6

Base case 0.18

402.8

−0.25

402.7

−1.35

Fin vertical position, Y2 (m) Y2 = −0.03 Y2 = 0.0 Y2 = −0.06 Y2 = −0.10

386.5 386.6 384.3 384.6

Base case 0.344 −1.76 −1.43

Fin spacing, S (m) S = 0.04 S = 0.035 S = 0.033

386.5 382.3 379

Base case −3.85 −7.32

Fin height, H (m) H = 0.7 H = 0.8 H = 0.9

386.5 377.2 373.5

Base case −9.13 −13.02

Active part vertical position, Y1 (m) Y1 = −0.17 386.5 Y1 = -0.10 387.2

Base case 1.26

Oil channel width, W (m) Base 0.001 m decrease in each channel 0.001 m increase in each channel 0.002 m increase in each channel

3. Optimization procedure In this section, the Response Surface Method will be utilized to optimize the geometry and find the minimum hotspot temperature. In this regard, the thermal behavior of the studied distribution transformer has been simulated according to the procedure described in the previous sections. Then, the temperature distribution, and thereby the value of hotspot temperature have been determined as the critical parameter for the optimization study. After that, the effects of various geometrical parameters on the hotspot temperature have been examined and the influential parameters have been identified. Using a specific procedure, which will be discussed later, the influential parameters are selected as the fin height, length and spacing for finding the optimized geometry with the minimum hotspot temperature. 3.1. Application of response surface method The response surface methodology (RSM) is a proficient mathematical and statistical technique useful for developing, optimizing, and improving products and industrial processes, allowing the determination of optimum parameters with a minimum number of runs. The RSM explores the relationships between the several explanatory variables and response variables, which was introduced by G. E. P. Box and K. B. Wilson in 1950s [35]. This method properly functions in designing the experiment and runs by defining the variables and also, the maximum and minimum limits of each variable. Response surface optimization is more advantageous than the traditional single parameter testing approaches that observe the effects of one factor at one time on the variable [36]. Additionally, RSM can be employed to develop empirical models that correlate the response to the effective factors. The most common types of RSM models are the three-level factorial, the face-centered central composite design (CCD) and the Box Behnken design. These symmetrical designs differ from each other with respect to their selection of experimental points, number of levels for variables, and number of runs and blocks. The CCD model was used in this study as it provides high-quality predictions over the entire design range [37]. Moreover, for finding optimum values of geometric parameters, the quadratic model of response surface methodology was utilized [38]. The required tests were identified using Design Expert, which is a powerful software for statistical simulation and analysis.

produce 282 cases for each the value of the hotspot temperature should be determined using ANSYS-Fluent software. Due to the complexity of the problem, each run in ANSYS-Fluent software takes nearly six days to converge; therefore, taking all parameters into account in RSM is so time consuming and computationally expensive. Therefore, initially, the most important parameters should be determined to reduce the number of required parameters. Thus, firstly, a parametric study was performed to determine the most effective parameters among these seven parameters. In this regards, considering the base case geometrical specifications, several values for each parameter were assigned, and for each set of the geometrical specifications with the help of numerical simulation, the hotspot temperature was determined. Therefore, the total number of 17 simulations were performed for this parametric study as reported in Table 7. It should be mentioned that to investigate the effect of a certain parameter on the hotspot temperature, all other parameters were set to the base case specifications. It is clear from Table 7 that three parameters including fin height, length and spacing can be chosen as the more influential parameters on the hotspot temperature. Findings these three parameters, the following ranges of the parameters were selected based on the geometry of the actual transformer. - Fin height: (0.5 m ≤ H ≤ 0.9 m) - Fin length: (0.06 m ≤ L ≤ 0.1 m) - Fin spacing: (0.032 m ≤ S ≤ 0.048 m)

3.2. Test method According to the configuration of the transformer illustrated in Figs. 1–4, various geometrical specifications may influence the hotspot of the transformer. Therefore, the main parameters can be listed as active part vertical position, oil channel width, fin width, fin height, fin length, fin vertical position and fin spacing. Considering all of these seven parameters, RSM program will

Moreover, five levels were considered for these parameters according to Table 8, and then the CCD model of RSM randomized and generated 20 run cases according to the number of parameters and their levels. 450

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Table 9 to generate the total of 20 cases by RSM. The hotspot temperature associated with each case is determined by running the numerical simulation. Table 9 indicates the details of all runs and the resulting hotspot temperature. In order to investigate the effects of these run conditions on the output variable (HST) and to find the sufficiency of the predicted model by RSM, the variance analysis (ANOVA) was implemented. This is a significant statistical tool to confirm the importance and accuracy of the generated models. Table 10 summarizes the ANOVA results, where for instance HL shows multiple of fin height and fin length and H2 refers to the fin height square. In this analysis, the interactions between the input variables have also been considered. As can be seen in Table 10, the statistical analysis of the model presents the value of 0.9904 for coefficient of multiple determinations (R2). As R2 is very close to 1 it states that this model is suitable to calculate the HST. However, R2 is not fully able to specify the accuracy of the model as it just shows the variations around the mean response; therefore, another coefficient, called the adjusted coefficient of determination (R2-Adj) was also used in this statistical analysis. In the calculation of R2- Adj, the sum of squares is substituted by the mean sum of squares which is more meaningful. In this model the value of 0.9817 was obtained for R2- Adj, which is less than R2; however, it is still very close to 1 and suggests an acceptable goodness-of-fit. F-value and P-value are two other significant factors that determine the accuracy and adequacy of the model. F-value is the ratio of the mean square of the model to the residual error. Calculated F-value should be greater than 4 if the model is a good predictor of the results. According to Table 10, the F-value of the present statistical model is 114.16 indicating that most of the variations in the response can be explained by the regression equation and the model is meaningful. On the other hand, the P-value is used to estimate whether F distribution is

Table 8 Design matrix in the CCD for the present problem. Parameter

Level

Fin height, H(m) Fin length, L(m) Fin spacing, S(mm)

-alpha

Low

Center

High

+alpha

0.5 0.06 32

0.6 0.07 36

0.7 0.08 40

0.8 0.09 44

0.9 0.1 48

As mentioned before, the response factor (dependent variable) in this model is the hotspot temperature (HST) of the transformer, and RSM method was used to find the optimum status of the input variables (independent parameters) to obtain the minimum value of HST. As mentioned before, RSM can also propose a correlation for the response factor as a function of the input variables and interaction effects of these variables. The multivariate model can be written as follows [39] K

Y = β0 +

∑ i=1

K

βi x i +

K

∑∑ i=1

K

βij x i x j +

j=1

∑ i=1

βii x i2 + ε

(10)

where, K is the number of independent variables, β0 denotes the constant coefficient, βi, βii and βij represent the linear coefficient of the ith factor, the quadratic coefficient of the ith factor, and the interaction of ith and jth factors, respectively. The ɛ refers to the residual associated with the experiments and Y is the dependent variable. The purpose is to find a correct relationship between the independent variables and the response. 4. Result and discussion 4.1. Temperature distribution Since, this study aims at presenting an optimized cooling system for the transformer through reducing the temperature of the hotspot, the details of temperature field and flow patterns, as well as the heat transfer mechanism in the transformer will not be discussed here. However, to gain a global view of the thermal distribution inside the transformer and in particular to become more familiarized with the possible locations of the hotspot, the oil temperature field in the central plane of the transformer Z = 0 (shown in Fig. 7a) is illustrated in Fig. 7b. It should be mentioned that there is also no fin in this plane as indicated in Fig. 7a. Fig. 7b demonstrated that the maximum temperature area is located on the top section of the narrow oil channels, where the main damage of the solid insulation of transformers is commonly reported. This observation is also consistent with the findings reported by other researchers [17,23,40]. As aforementioned, the main role of the oil channels is to transfer the generated heat in the active parts to the tank surfaces and fins by natural convection, where the heat is rejected to the surrounding area. Clearly, the oil channels are of great importance in this process, however, their performances are limited by the relatively narrow flow area and the viscosity of the oil, which cause the oil temperature reaches its maximum value in these channels. Fig. 7b also shows that temperature in the upper section is higher than the lower section, which is due to the buoyancy force in the natural convection. 4.2. Statistical analysis (ANOVA) In this section, the relationship between the independent variables and the response that is the hotspot temperature is estimated using the design of experiment method, according to the RSM. On this basis, three main design parameters of fin height, length and spacing are considered as the independent variables as discussed earlier. Different combinations of these variables were also identified as presented in

Fig. 7. a) Schematic of Z = 0 plane in transformer; b) Temperature contour and HST of transformer on Z = 0 plane. 451

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

general form of the relationship between the effective parameters (height, length and spacing of fin) and the response (hotspot temperature) can be presented by Eq. (11). Considering Table 10 the effects of all linear terms with the P-value of < 0.05 are highly significant. The quadratic terms of H and L are not significant, as their P-values fall above 0.05. However, the significance of S2 with the P-value lower than 0.05 is found to be higher than the other quadratic factors. Additionally, all the interaction factors have P-values higher than 0.05 have no practical significance. It is worth noting that Eq. (11) is a very valuable and useful correlation that is capable of predicting the hotspot temperature based on three most significant geometrical parameters of the cooling system.

Table 9 Values of independent variables generated by RSM along with their responses. Run

H(m)

L(m)

S(m)

Response: HST (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

0.6 0.8 0.7 0.7 0.8 0.6 0.6 0.7 0.8 0.8 0.9 0.7 0.6 0.5 0.7 0.7 0.7 0.7 0.7 0.7

0.09 0.07 0.08 0.08 0.09 0.09 0.07 0.08 0.09 0.07 0.08 0.08 0.07 0.08 0.08 0.08 0.08 0.08 0.06 0.1

0.036 0.044 0.040 0.032 0.036 0.044 0.036 0.040 0.044 0.036 0.040 0.040 0.044 0.040 0.040 0.040 0.040 0.048 0.040 0.040

381.49 387.33 386.40 378.61 370.67 390.61 392.47 386.40 377.37 379.60 372.81 386.40 402.85 401.25 386.40 386.40 386.40 389.57 397.96 376.55

HST = 386.23−6.83H −5.31L + 3.49S−0.66S 2

4.5. Perturbation analysis Another useful statistical analysis is to show the effects of all factors on the response by the perturbation plot, as shown in Fig. 9. The perturbation plot shows how the response changes as each factor moves from the chosen reference point, with all other factors held constant at the reference value. Design-Expert sets the reference point default at the middle of the design space (the coded zero level of each factor). In this plot the slope of each parameter line shows the sensitivity of the response to that parameter. Clearly, the HST is strongly influenced by the fin height as indicated by the slop of the A line. Furthermore, the increase of factors of A (fin height) and B (fin length) lead to the reduction of HST, which is favorable, whereas, the hotspot temperature increases with the increase in fin spacing as denoted by the factor C in Fig. 9. Fig. 10 displays the 3D plots of HST versus a) H and L, b) H and S, and c) L and S. Since, this kind of graph shows the interaction of two independent parameters in their whole range in a single graph, it is very meaningful and practical. Fig. 10a clearly indicates that in the whole range of L and H, increase of one parameter while the other one is constant, leads to diminishing of HST. Additionally, from Fig. 10b and 10c, it is observed that by increasing the value of S in the constant values of H or L, HST undergoes a raise that is unfavorable in cooling performance of transformers.

Table 10 The variance analysis (ANOVA) results of the HST. Source

Sum of Square

Model 1412.09 H 746.83 L 450.55 S 194.83 HL 2.33 HS 3.23 LS 0.66 2 H 0.15 L2 0.45 S2 10.89 Residual 13.74 R2:0.9904 R2- Adj:0.9817

(11)

df

Mean Square

F-value

P-value

9 1 1 1 1 1 1 1 1 1 10

156.90 746.83 450.55 194.83 2.33 3.23 0.66 0.15 0.45 10.89 1.37

114.16 543.40 327.82 141.76 1.70 2.35 0.48 0.11 0.32 7.92

< 0.0001 < 0.0001 < 0.0001 < 0.0001 0.2219 0.1562 0.5046 0.7487 0.5817 0.0183

statistically significant or not, and the model is statistically significant, when the P-value is lower than 0.05. Since the P-value of the present model is < 0.001, it can be concluded that the model is statistically significant. 4.3. Diagnostic analysis

4.6. Optimization of response

In this section, some diagnostic plots are presented to investigate the statistical qualification of the model. The most important diagnostic plot is the normal probability plot of the studentized residuals. The data points in this plot should be approximately linear; a non-linear pattern indicates non-normality in the error term. As Fig. 8a indicates, there is no sign of problem in data, and points on normal probability plot conform to a straight line. Predicted values versus actual values plot are displayed in Fig. 8b, which also indicate a good correlation between the predicted and actual data. Moreover, Fig. 8c illustrates the residuals of runs, which are in an acceptable range.

One of the most important goals of response surface method is finding an optimum response [5]. In searching for an optimum response, the desirability, ranged from zero to one, plays a predominant role to show whether the predicted response is optimized or not. Value of one represents the ideal optimized response, while zero indicates that the predicted response has a poor quality and is not an optimum response. To this end, the response surface optimization is an ideal technique for determination of the best geometry combination. Here, the goal is to minimize the hotspot temperature of the transformer to reduce the risk of failures and defects. The optimum geometrical parameters obtained by RSM are the fin height of 0.8 m, fin length of 0.09 m and fin spacing of 0.036 m, as shown in Table 12. Considering this optimum geometry, the predicted minimum hotspot temperature is 370.8 K, which is about 16 °C lower than the hotspot temperature of the actual studied transformer. This reduction is very significant, since it has been reported that the rate of the Loss-of-Life of transformer exponentially increases with the temperature raise. Based on a semi-empirical correlation reported in [41], the expected life of the optimized distribution transformer that runs with the hotspot temperature of 370.8 K is about 6 times the actual transformer, which runs with the hotspot temperature of 386.4 K.

4.4. Model estimation The RSM is capable of proposing a correlation for the relationship between the response factor and input variables. According to Table 11, a combination of different relations between the response and input variables can be established such as a linear relation, quadratic relation and interactions of variables. The constant coefficients of all these possible terms as estimated by RSM are listed in Table 11. To determine the significance of each term with respect to other terms their corresponding P-values from Table 10 must be considered. Accordingly, the 452

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Fig. 8. a) Predicted values versus actual values; b) Normal Probability; c) Residual versus run. Table 11 Estimated coefficients for HST. Factor

Coefficient Estimate

Intercept H L S HL HS LS H2 L2 S2

386.233 −6.8320 −5.306 3.489 0.539 −0.635 −0.286 0.077 0.133 −0.658

5. Conclusion Three dimensional thermal analysis of a 200kVA distribution transformer has been performed numerically using a CFD software in search for an optimized geometrical cooling system. The thermal variations of the thermo-physical properties of the transformer oil are determined experimentally and incorporated in the numerical modeling. Response surface method has been employed as the optimization tool, which requires the effective geometrical parameters as the input variables and the transformer hotspot temperature as the response.

Fig. 9. Perturbation graph for HST.

Through numerical parametric study the most influential geometrical parameters on the thermal performance of the transformer are identified as the fin height, length and spacing. Application of the RSM 453

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

Fig. 10. Interaction effects of a) H and L b) H and S c) L and S on HST.

transformer, which ensures better performance and increases the life expectancy of the transformer greatly.

Table 12 Response optimization for hotspot temperature. Parameters

HTS

Goal

minimum

Optimum combination

Height 0.8 m

Length 0.09 m

Spacing 0.036 m

Predicted response

Desirability

370.8 K

1

Acknowledgements The authors gratefully acknowledge the Khorasan Regional Electricity Company for the financial support (contract No. 101280), and Dr. Abdollah Kamyab for providing useful information and advice on various technical issues.

optimization procedure dictates a number of runs with various combinations of the influential geometric parameters, through which the optimal geometrical design, which minimizes the transformer hotspot temperature are determined. Furthermore, RSM statistical analysis can provide information about how the hotspot temperature is affected by the geometrical specifications. In fact, a correlation for the variation of hotspot temperature as a function of fin height, length, and spacing has been proposed, which indicates the usefulness of the RSM optimization procedure in practical applications. Results indicate that hotspot temperature is more influenced by the fin height as compared to the fin length and spacing, yet their effects on the hotspot temperature are different. As expected, increasing in fin height and length reduces the hotspot temperature, while reducing the fin spacing, which provides more heat transfer area acts in a favorable manner on the hotspot temperature. The optimum geometry that minimizes the hotspot temperature is found to have the fin height of 0.8 m instead of 0.7 m, fin length of 0.09 m instead of 0.08 and fin spacing of 0.036 m instead of 0.04 m. The hotspot temperature of this optimized geometry is about 16 °C lower than that of the actual

References [1] Sefidgaran M, Mirzaie M, Ebrahimzadeh A. Reliability model of the power transformer with ONAF cooling. Elect Power Energy Syst 2012;35(1):97–104. [2] Schneider J, et al. Asset management techniques. Elect Power Energy Syst 2006;28(6):643–54. [3] Krause C. Power Transformer Insulation – History, Technology and Design. IEEE Trans Dielectr Electr Insul 2012;19(6):1941–7. [4] Perez, J, Fundamental Principles of Transformer Thermal Loading and Protection, in Protective Relay Engineers; 2010, IEEE: College Station, TX. [5] Fu W, McCalley JD, Vittal V. Risk assessment for transformer loading. IEEE Trans Power Syst 2001;16(3):346–53. [6] Taghikhani MA, Gholami A. Prediction of hottest spot temperature in power transformer windings with non-directed and directed oil-forced cooling. Elect Power Energy Syst 2009;31(7–8):356–64. [7] Pradhan MK, Ramu TS. Prediction of hottest spot temperature (HST) in power and station transformers. IEEE Trans Power Delivery 2003;18(4):1275–83. [8] Susa D, Nordmann H. A simple model for calculating transformer hot-spot temperature. IEEE Trans Power Del 2009;24(3):1257–65. [9] Hajidavalloo E, Mohamadianfard M. Effect of sun radiation on the thermal behavior of distribution transformer. Appl Therm Eng 2010;30(10):1133–9. [10] Souza LM, et al. Thermal modeling of power transformers using evolving fuzzy

454

Electrical Power and Energy Systems 104 (2019) 443–455

L. Raeisian et al.

[26] Tsilia MA, et al. Power transformer thermal analysis by using an advanced coupled 3D heat transfer and fluid flow FEM model. Int J Therm Sci 2012;53:188–201. [27] Chereches NC, et al. Numerical study of cooling solutions inside a power transformer. Energy Procedia 2017;112:314–21. [28] Eslamian M, Vahidi B, Eslamian A. Thermal analysis of cast-resin dry-type transformers. Energy Convers Manage 2011;52(7):2479–88. [29] Gastelurrutia J, et al. Numerical modelling of natural convection of oil inside distribution transformers. Appl Therm Eng 2011;31(4):493–505. [30] Cengel YA, Ghajar AJ, Heat and mass transfer: Fundamentals & applications. 5 ed; 2015, New YorK: McGraw-Hill Education. [31] Roohoulamini B, Heat Transfer Investigation of Air-type Electric Transformers and the Effect of Parameters on Their Performance Using a Computational Method in Topical Cities of Kerman Province, in Mechanical engineering department2016, MSc thesis, Kerman University of Advanced Technology: Kerman, Iran. [32] Bar-Cohen A, Rohsenow WM. Thermally optimum spacing of vertical, natural Convection Cooled, Parallel Plates. J Heat Transfer 1984;106(1):116–23. [33] Churchill SW, Usagi R. A general expression for the correlation of rates of transfer and other phenomena. AIChE J 1972;18(6):1121–8. [34] Versteeg, HK, Malalasekera W, Introduction to Computational Fluid Dynamics. 1st ed; 1995: Wiley. [35] Steinberg DM. George box and the design of experiments: statistics and discovery. Appl Stochastic Models Bus Ind 2014;30(1):36–45. [36] Awad OI, et al. Response surface methodology (RSM) based multi-objective optimization of fusel oil-gasoline blends at different water content in SI engine. Energy Convers Manage 2017;150:222–41. [37] Danmaliki GI, Saleh TA, Ahmad SA. Response surface methodology optimization of adsorptive desulfurization on nickel/activated carbon. Chem Eng J 2017;313(1):993–1003. [38] Avramovic, JM, et al., Optimization of Sunflower Oil Ethanolysis Catalyzed by Calcium Oxide: RSM Versus ANN-GA. Energy Conversion and Management; 2015. 105: p. 1149–1156. [39] Betiku E, et al. Predictive capability evaluation of RSM, ANFIS and ANN: a case of reduction of high free fatty acid of palm kernel oil via sterification process. Energy Convers Manage 2016;124:219–30. [40] Campos ART, Mariscal IC, Hernandez SG. Simulation of a distribution transformer. WSEAS Trans Fluid Mechan 2012;7(3):106–15. [41] Brown ER, Electric power distribution reliability. CRC Press: 2008,.

systems. Eng Appl Artif Intell 2012;25(5):980–8. [11] Koufakis EI, Karagiannopoulos CG, Bourkas PD. Thermal coefficient measurements of the insulation in distribution transformers of a 20 kv network. Measurement 2008;41:10–9. [12] Abu-Elanien AEB, Salama MMA, A Monte Carlo Approach for Calculating the Thermal Lifetime of Transformer Insulation. Elect Power Energy Syst; 2012 43 (1): p. 481–487. [13] Hosseini R, Nourolahi M, Gharehpetian GB. Determination of OD cooling system parameters based on thermal modeling of power transformer winding. Simul Model Pract Theory 2008;16:585–96. [14] Torriano F, Chaaban M, Picher P. Numerical study of parameters affecting the temperature distribution in a disc-type transformer winding. Appl Therm Eng 2010;30(14–15):2034–44. [15] Zhang J, Li X. Oil cooling for disk-type transformer windings—Part II: parametric studies of design parameters. IEEE Trans Power Deliv 2006;21(3):1326–32. [16] Torriano F, Picher P, Chaaban M. Numerical investigation of 3D flow and thermal effects in a disc-type transformer winding. Appl Therm Eng 2012;40:121–31. [17] Campelo, H., et al., Detailed CFD Analysis of ODAF Power Transformer, in International Colloquium Transformer Research and Asset Management 2009, cigre: Croatia. [18] Skillen A, et al. Numerical prediction of local hot-spot phenomena in transformer windings. Appl Therm Eng 2012;36:96–105. [19] Coddé J, Veken WV, Baelmans M, Assessment of a Hydraulic Network Model for Zigezag Cooled Power Transformer Windings. Appl Therm Eng 80: p. 220–228. [20] Yatsevsky VA. Hydrodynamics and heat transfer in cooling channels of oil-filled power transformers with multicoil windings. Appl Therm Eng 2014;63:347–53. [21] Rahimpour E, Barati M, Schäfer M. An investigation of parameters affecting the temperature rise in windings with zigzag cooling flow path. Appl Therm Eng 2007;27(11–12):1923–30. [22] Smolka J. CFD-based 3-D optimization of the mutual coil configuration for the effective cooling of an electrical transformer. Appl Therm Eng 2013;50(1):124–33. [23] El Wakil N, Chereches N-C, Padet J. Numerical study of heat transfer and fluid flow in a power transformer. Int J Therm Sci 2006;45(6):615–26. [24] Oh K-J, Ha S-S. Numerical calculation of turbulent natural convection in a cylindrical transformer enclosure. Asian Res 1999;28:429–41. [25] Radakovic Z, Jevtic M, Das B. Dynamic thermal model of kiosk oil immersed transformers based on the thermal buoyancy driven air flow. Elect Power Energy Syst 2017;92:14–24.

455