Prediction of hottest spot temperature in power transformer windings with non-directed and directed oil-forced cooling

Electrical Power and Energy Systems 31 (2009) 356–364

Contents lists available at ScienceDirect

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

Prediction of hottest spot temperature in power transformer windings with non-directed and directed oil-forced cooling M.A. Taghikhani *, A. Gholami Electrical Engineering Department, Iran University of Science and Technology (IUST), Narmak, 16846 Tehran, Iran

a r t i c l e

i n f o

Article history: Received 30 January 2008 Received in revised form 24 February 2009 Accepted 16 March 2009

Keywords: Power transformer Temperature distribution Hottest spot Forced cooling Heat equation

a b s t r a c t Power transformer outages have a considerable economic impact on the operation of an electrical network. One of the most important parameters governing transformer’s life expectancy is the hottest spot temperature (HST) value. The classical approach has been to consider the hottest spot temperature as the sum of the ambient temperature, the top-oil temperature rise, and the hottest spot to top-oil temperature gradient. The authors proposed a numerical method based on heat transfer theory using the finite element method and they only needed to solve heat conduction equation. The transformer selected for simulation was a 32 MVA transformer with non-directed oil-forced (NDOF) cooling and directed oil-forced (DOF) cooling. A comparison of the authors results with those obtained from finite integral transform and experimental test confirms the validity and accuracy of the proposed method. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Large power transformers are amongst the most valuable assets in electrical power networks and an outage affects the stability of the network. In the power transformer a part of the electrical energy is converted to heat. Although this part is quite small compared to the total electric power transferred through a transformer, it causes significant temperature rise. The thermal impact leads not only to long-term oil/paper-insulation degradation; it is also a limiting factor for the transformer operation [1]. Therefore, the knowledge of temperature, especially the hottest spot temperature, is of high interest. The basic condition for transformer loading is the hottest spot temperature. In order to avoid insulation failures it must not exceed the prescribed value. The hottest spot calculation procedure is given in the International Standards [2–4]. An algorithm has been given for calculating the hottest spot temperature of a directly loaded transformer in [5,6]. These papers propose improvements in the modeling of thermal processes inside the transformer tank. A simple equivalent circuit to represent the thermal heat flow equations for power transformers is presented in [7]. The use of an average heat transfer coefficient is typical in the transformer design process used to calculate the required number (area) of

* Corresponding author. Tel.: +98 21 77240492/77240493x2600; fax: +98 21 77240490. E-mail addresses: [email protected] (M.A. Taghikhani), [email protected] (A. Gholami). 0142-0615/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2009.03.009

cooling surfaces [8–13]. An electromagnetic analysis using a finite element model is used to calculate the eddy losses in the transformer windings and only equivalent circuit to represent the thermal model has been established in [14]. Finite element analysis with ANSYS software to calculate the transformer losses and temperature distribution in the windings is proposed in [15], but oil temperature around the windings is considered constant and is modeled only layer windings. Fluid flow analysis with FEMLAB software is used in [16,17] to calculate temperature distribution in the transformer and cooling system but this method is time consuming and impossible in some cases. The authors note that the IEEE loading guide and other similar documents offer relations for the calculation of the hottest spot temperature based on per-unit load. The formulations tend to ignore the possibility of two transformers that have identical ratings but have a different winding structure and varying heat loss/unit volume. The method suggested by the authors gives due representation for this difference and, hence, is believed to give more accurate estimates. In this paper, a procedure is proposed for obtaining the temperature distribution in the power transformer windings. The procedure requires calculation of the heat transfer coefficients in the heat conduction equation. The finite element method is used in the numerical solution. For this purpose, we have provided a code using MATLAB version 6.5. In this paper, the heat conduction theory is formulated in Section 2 while its solution using finite element method is presented in Section 3. Conclusions are provided in Section 4.

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

357

Nomenclature a b Cp Dm f1 f2 f3 f4 fr g Gr Grhf Gz h1 h2 h3 h4 K k1 k2 k3 k4 kcu kin koil kr kz l m n Nu Num Pr Q Q0 qw

inner radius of the annular disc or layer (m) outer radius of the annular disc or layer (m) specific heat of transformer oil (J kg1 K1) mean diameter of annular disc or layer of windings (m) boundary function at inner cylinder surface boundary function at outer cylinder surface boundary function at bottom cylinder surface boundary function at top cylinder surface relative index of free to forced convection acceleration due to gravity (m s2) Grashof number based on temperature difference Grashof number based on constant heat flux Graetz number heat transfer coefficient at inner surface (W m2 K1) heat transfer coefficient at outer surface (W m2 K1) heat transfer coefficient at bottom surface (W m2 K1) heat transfer coefficient at top surface (W m2 K1) referred to as equivalent thermal conductivity thermal conductivity at inner cylinder surface (W m1 K1) thermal conductivity at outer cylinder surface (W m1 K1) thermal conductivity at bottom cylinder surface (W m1 K1) thermal conductivity at top cylinder surface (W m1 K1) thermal conductivity of copper (W m1 K1) thermal conductivity of insulation (W m1 K1) thermal conductivity of transformer oil (W m1 K1) thermal conductivity in radial direction (kr = k1 = k2) thermal conductivity in axial direction (kz = k3 = k4) thickness of the disc or height of layer winding (m) oil temperature gradient along winding height (K m1) temperature rise exponent to bottom oil local Nusselt number mean Nusselt number Prandtl number volumetric heat source function (W m3) volumetric heat source at ambient temperature (W m3) heat flux in (W m2)

2. Mathematical formulation for heat conduction equation The geometry of the transformer winding is complex and is not typical. Under general conditions, we assume that the transformer winding is cylindrical; hence, a layer or a disc winding is a finite annular cylinder [10]. The thermal and physical properties of the system would be equivalent to a composite system of insulation and conductor. It has been assumed that heat is generated throughout body at a constant rate, and oil in the vertical and horizontal ducts take away the heat through the process of convection. In an actual transformer winding, the conductor is the only heat source. Later in this section, formulations are given for calculating different thermal and physical properties of the system. It has been assumed that temperature is independent of space variable, because the winding structure is symmetrical. The surface of the disc or layer has been assumed to be flat. The generalized system of non-homogeneous heat conduction equations with non-homogeneous boundary conditions, in Cartesian coordinate system are written [10,18]

R Ra Rahf Re ri’s Tamb Tb tcu tin Tm Ttop

loss ratio = load loss/no-load loss Rayleigh number based on temperature difference Rayleigh number based on constant heat flux Reynolds number i = 1,2,. . .radius of insulation and conductor layers (m) ambient temperature (K) temperature at the bottom of the winding (K) thickness of conductor (m) thickness of insulation (m) temperature average of oil and winding surface (K) temperature at the top of the winding (K)

Abbreviations DOF directed oil-forced cooling HST hottest spot temperature HTC heat transfer coefficient HV high voltage LV low voltage NDOF non-directed oil-forced cooling p.u. per unit Greek symbols b coefficient of volumetric expansion of oil (K1) ultimate bottom oil temperature rise over ambient temhu perature(K) ultimate bottom oil temperature rise over ambient temhfl perature at rated load (K) d characteristic dimension (m) q density of transformer oil (kg m3) f(n) correction factor to cylindrical curvature ac temperature coefficient of electrical resistance (K1) l viscosity of oil (kg m1 s1) lb viscosity of oil at oil bulk mean temperature (kg m1 s1) lw viscosity of oil at winding wall temperature (kg m1 s1) Subscripts b bulk property m mean value w wall value

@2T @2T k þ @x2 @y2

! þQ ¼0

ð1Þ

in the region a < x < b and 0 < y < l. At the inner cylindrical surface

k1

@T þ h1 T ¼ f1 ðyÞ @x

ð2Þ

At the outer cylindrical surface

k2

@T þ h2 T ¼ f2 ðyÞ @x

ð3Þ

At the bottom flat surface

k3

@T þ h3 T ¼ f3 ðxÞ @y

ð4Þ

At the top flat surface

k4

@T þ h4 T ¼ f4 ðxÞ @y

ð5Þ

358

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

Eqs. (1)–(5) represent the general heat conduction equation with convection at all four boundary surfaces. In the above equations, temperature T is a function of space variables x and y. The term Q is the heat source function and has been modified here to take care of variation of electrical resistance of copper with temperature. The heat source term Q is

Q ¼ Q 0  ½1 þ ac  ðT  T amb Þ

ð6Þ

across the four surfaces. It is reported [10], that the heat transfer coefficient depends on as many as 13 factors (e.g. winding size, type, duct dimensions, oil velocity, type of oil circulation, heat flux distribution, oil thermal properties, etc.). In [10], corrections have been given for temperature dependence of the thermal and physical properties of oil, such as viscosity, specific heat, volumetric expansion, and thermal conductivity. We used some of the heat transfer relations in natural cooling (ON) mode to calculate the heat transfer coefficients [18]. The local Nusselt and Rayleigh numbers for laminar flow over vertical plates are

ac is the temperature coefficient of electrical resistance of the copper wire. Therefore, the function Q becomes a temperature dependent, distributed heat source. Boundary functions f1(y) and f2(y), derived from Newton’s law of cooling are of the form (7) and (8)

Nu ¼ 0:6Ra0:2 hf

ð12Þ

f1 ðyÞ ¼ h1  ðT b þ m1 :yÞ

ð7Þ

Rahf ¼ Grhf  Pr

ð13Þ

f2 ðyÞ ¼ h2  ðT b þ m2 :yÞ

ð8Þ

The term Tb is the temperature at the bottom of the disc or layer, as applicable. Terms m1 and m2 are the temperature gradient along the winding height (for layer) or along the disc thickness for a disc. Similarly, functions f3 and f4 represent temperature profiles across bottom and top surfaces, and have the same form as shown in (7) and (8), where temperature gradient term has been taken as zero.

f3 ðxÞ ¼ h3  T b

ð9Þ

f4 ðxÞ ¼ h4  T top

ð10Þ

Fig. 1 shows boundary conditions of a winding. Thermal conductivities are unequal in different directions. But in an actual case, the thermal conductivities in radial directions (k1 and k2) are equal and conductivities in the axial direction will also be the same (k3 and k4). Thermal conductivity has been treated as a vector quantity. The resultant thermal conductivity of the system can be estimated as

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 K ¼ kr þ kz

ð11Þ

where Rahf and Grhf are the local Rayleigh and Grashof numbers based on heat flux (qw) at characteristic dimension (d). Pr is the Prandtl number. The Rayleigh number based on constant heat flux is

Rahf ¼

kr ¼ 

r logr2 1

k1

kz ¼

þ

r logr3 2

Num ¼ 1:5½Nud¼l

k2

þ ::: þ

logr rn

n1

kn

kcu  kin  ðt cu þ tin Þ tin  kcu þ t cu  kin

Term K represents the resultant thermal conductivity of the insulation and conductor system. Heat transfer coefficients (HTC) h1–h4, are different across all four surfaces. To determine boundary functions f1–f4, it is necessary to calculate the heat transfer coefficients

ð15Þ

However, the correction factor to be applied to formula (12) for cylindrical curvature is (30 < Pr < 50):

f ðnÞ  1 þ 0:12  n

ð16Þ

where

n¼

pffiffiffi 2 2 Gr0:25



d r

ð17Þ

Here Gr is the Grashof number based on the temperature difference. The heat transfer coefficient can be computed as

h¼ 

ð14Þ

2

koil  l

The mean Nusselt number in this case can be computed as

where

log rrn1

g  b  C p  q2  qw  d4

Nu  koil d

ð18Þ

where h represents local coefficient. The mean coefficient (hm) can be calculated from the mean Nusselt number, as in (15). Knowing hm for a particular surface, the temperature difference between the winding surface and oil can be found by dividing hm by the heat flux through the surface. The mean Nusselt number of the top surface of the annular cylindrical winding for the laminar and turbulent regime is

Num ¼ 0:54 Ra0:25 ¼ 0:61 Ra0:2 hf

ð19Þ

1 4

1 3

Num ¼ 0:15 Ra ¼ 0:24 Rahf d¼

ð20Þ

ba 2

ð21Þ

where a and b are the inner and outer radius of the annular disc or the layer. The Nusselt number of the bottom surface of the annular cylindrical winding for the laminar and turbulent regime is given in (22)

Num ¼ 0:27 Ra0:25 ¼ 0:35 Ra0:2 hf

ð22Þ

The axial oil temperature gradient in the presence of cooling by a constant heat flux can be found by using Eq. (23), from [10] thus

@T q ¼ 4:25  102 w  @y koil

Fig. 1. Boundary conditions for a winding.



l p  Dm

49

1

Rahf9

ð23Þ

where l is the winding height, Dm is the mean diameter of the annular disc of the windings. Determination of boundary conditions in forced convection (OF mode) too requires calculating HTC. The

359

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

h i1  l 0:14 1 4 3 b Num ¼ 1:75  Gz þ 0:012  ðGz  Gr3 Þ3 

lw

h i 1 Num ¼ 1:63  1:8  Gz0:9 þ 0:02  ðGz  Gr3hf Þ1:16  1 3

D l

Gz ¼ Re  Pr 

ð24Þ



lb lw

0:125 ð25Þ ð26Þ

Gz is called Graetz number. The terms lb and lw are viscosity of oil computed at oil bulk mean temperature of oil and at winding wall temperature, respectively. Re is the Reynolds number. The relative importance of natural and forced cooling is indicated by the factor fr = Gr/Re2. If fr  1 then both of the cooling modes have to be considered. At a lower value of this factor, natural cooling can be ignored. The oil viscosity is an important property, which depends on temperature, is expressed in (27) from [9,10]

l ¼ a  exp



c

 ð27Þ

Tm

where

a = 0.0000013573 (kg m1 s1) c = 2797.3 (K) Fig. 2. Geometrical shape of the modeled windings.

expression when the oil velocity is lower values, the mean Nusselt number based on temperature difference is of the form of (24), corresponding mean Nusselt number based on constant heat flux is of the form of (25)

The viscosity was calculated at the mean oil and wall temperature. Initially, the winding surface temperature is not known, so a starting guess for the winding surface temperature needs to be made, after calculating the value of h, the temperature difference (Tw–Toil) must agree with the assumed value. To calculate the winding wall temperature at different surfaces, only the bottom oil temperature is necessary. In this paper, the bottom oil rise above ambient temperature has been calculated as

Table 1 Characteristics of transformer windings. Transformer rating (MVA)

Winding type

Cu-loss (W)

Dimensions (mm) a,b,l

Remarks

32 32 32

LV layer1 LV layer2 HV disc

7579 8536 570

263,287.5,1460 301.5326,1460 371,438,21

Per layer Per layer Per disc

122 120 Temperature (°C)

118 116 114 112 110 108 106

1

0.285 0.28 0.275

0.5 0.27 Height of LV layer1 winding(m)

0

0.265 Radial thickness(m)

Fig. 3. Temperature distribution of LV layer1 winding at Tamb = 25 (°C) with Re = 750.

360

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

Temperature(ºC)

120

115

110

105

0.325

1 0.32 0.315

0.5

Height of LV layer2 winding(m)

0.31 0

0.305 Radial thickness(m)

Fig. 4. Temperature distribution of LV layer2 winding at Tamb = 25 (°C) with Re = 750.

!n I2r  R þ 1 hu ¼ hfl  Rþ1

ð28Þ

where hu is the bottom oil temperature rise over ambient temperature, hfl is the full load bottom oil temperature rise over ambient temperature obtained from an off-line test, and R is the ratio of load loss at rated load to no-load loss. The variable Ir is ratio of the specified load to rated load:

Ir ¼

I Irated

ð29Þ

Fig. 5. NDOF cooling in the HV discs.

The exponent n depends upon the cooling state. The IEEE loading guide recommends the use of n = 0.8 for natural convection and n = 0.9–1.0 for forced cooling. Assumptions have been made in the calculation of HTC in the ducts provided in the disc-type windings under oil-forced (OF) modes of the heat transfer (DOF and NDOF). The cooling in OF mode is due to mixed mode (natural and forced) convection. While the oil-flow velocity in both vertical and horizontal ducts has been assumed equal, we are in the DOF mode. If the flow velocity in horizontal ducts is assumed negligible compared to velocity in the vertical ducts we are in NDOF mode. In case of DOF mode, the HTC was assumed a function of both heat flux through the surface and the oil-flow velocity. The mechanism of heat transfer in this mode of

Fig. 6. DOF cooling in the HV discs.

361

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

cooling is same in both axial and radial direction of the disc. In case of NDOF mode, the HTC in the vertical duct has been estimated by the same formula as for DOF. However, the convection of heat in the horizontal ducts has been assumed a purely natural type.

3. Results and discussion We have provided a code using MATLAB version 6.5 for the finite element solution of the heat conduction equation. The model

110

Temperature(ºC)

108 106 104 102 100 98 0.02 0.015

0.43 0.42

0.01

0.41 0.4

0.005 Height of HV disc winding(m)

0.39 0.38

0

Radial thickness(m)

Fig. 7. Temperature distribution of HV disc winding in NDOF mode at Tamb = 25 (°C) with Re = 750.

82 81.5

Temperature(ºC)

81 80.5 80 79.5 79 78.5 78 0.02 0.015

0.43 0.42

0.01

0.41 0.4

0.005 Height of HV disc winding(m)

0.39 0.38

0

Radial thickness(m)

Fig. 8. Temperature distribution of HV disc winding in DOF mode at Tamb = 25 (°C) with Re = 750.

Table 2 HST (°C) magnitudes at Tamb = 25 (°C) with Re = 750. Load (p.u.)

0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

HST (°C) – LV layer1

95 108 123 140 158 177 199 222

HST (°C) – LV layer2

94 107 122 138 156 175 197 220

HST (°C) – HV disc NDOF

DOF

87 99 112 126 141 157 174 192

67 74 82 91 101 112 124 137

Global HST (proposed)

HST-Ref. [10]

95 108 123 140 158 177 199 222

92 104 117 134 152 183 200 229

362

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

has been validated on a transformer of rating 32 MVA. Fig. 2 shows the geometrical shape of the modeled windings. The characteristics of the above transformer windings have been tabulated in Table 1. The result of temperature distribution for LV layer1 winding in 1 p.u. load with Re = 750 has been shown in Fig. 3. It can be pointed out that for LV layer1 winding, maximum temperature location is around 90–95% of winding height from the bottom and at about 50% of radial thickness of the layer. The result of temperature distribution from LV layer2 winding has been shown in Fig. 4. It can be pointed out that for LV layer2 winding, maximum temperature location is the same as LV layer1 winding. Figs. 5 and 6 show NDOF and DOF cooling in the HV discs. Temperature distribution from a disc of HV winding has been shown in Figs. 7 and 8. It may be observed that the maximum temperature occurs in the neighborhood of 30–35% of the axial and 50% of the radial thickness of the disc in NDOF mode and 50–55% of the axial

and 50% of the radial thickness of the disc in DOF mode. Magnitudes of HST at different loading and OF modes are given in Table 2. In this work, different rates of forced oil circulation are taken into account. Figs. 9–12 show variations of local HST with load and Reynolds number in different windings and different OF modes. Ref. [10] has proposed an analytical method based on the boundary value problem of heat conduction in a transformer winding using finite integral transform techniques. This technique uses Fourier transforms and Hankel transforms and calculates the temperature based on polynomials with infinite order (sequential series) of space variables and time, and then approximates the temperature with finite order. The authors in [10] consider that by increasing the order of the polynomial for accuracy, the results will be more accurate, but it is not true. On the other hand, if we increase the order of the polynomial, the solution time will rise

260 240 220

180 160 140 120 100

1.4 1.3

1400 1200

1.2 1000

1.1

800

1

Load(p.u.)

600 0.9

400 0.8

Reynolds number

200

Fig. 9. HST (local) distribution of LV layer1 winding at Tamb = 25 (°C).

260 240 220 200 HS T(ºC)

HST(ºC)

200

180 160 140 120 100

1.4 1.3

1400 1200

1.2 1000

1.1

800

1 Load(p.u.)

600 0.9

400 0.8

200

Rey nolds num ber

Fig. 10. HST (local) distribution of LV layer2 winding at Tamb = 25 (°C).

363

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

180

HS T(ºC)

160

140

120

100

1.4 1.3

1400 1200

1.2 1000

1.1

800

1

Load(p.u.)

600 0.9

400 0.8

Rey nolds num ber

200

Fig. 11. HST (local) distribution of HV disc winding in NDOF mode at Tamb = 25 (°C).

140 130

HS T(ºC)

120 110 100 90 80 70

1.4 1.3

1400 1200

1.2 1000

1.1 Load(p.u.)

800

1

600 0.9

400 0.8

200

Rey nolds num ber

Fig. 12. HST (local) distribution of HV disc winding in DOF mode at Tamb = 25 (°C).

and the analytical method will not differ from numerical method. The analytical method in [10] is suitable for cylindrical windings with rectangular cross-section, but this limitation is not applicable for the finite element method. Tables 2 and 3 show comparison of the proposed method with boundary value method using finite integral transform used in [10] and experimental test in [9]. It is clear from Table 2 that the HST in

the case of NDOF mode is greater than that of corresponding DOF mode so is the average temperature rise. From Fig. 12 it is also borne out that by an increase in the amount of cooling by way of increased oil flow speed would provide a reduced winding temperature, which is obvious. However, an increase oil flow has a lesser effect on HST in case of NDOF mode than the DOF mode (Figs. 11 and 12).

364

M.A. Taghikhani, A. Gholami / Electrical Power and Energy Systems 31 (2009) 356–364

Table 3 HST (°C) locations (in LV winding) at Tamb = 25 (°C). Load (p.u.)

Proposed (radial thickness % – height %)

Ref. [9] (radial thickness % – height %)

0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

50/95 50/94.8 50/94.6 50/94.5 50/94.5 50/94.3 50/94.2 50/94

44/94 45/95 45/94 45/93 45/93 49/91 49/91 52/90

4. Conclusions In this paper, an attempt has been made to suggest a method to improve the accuracy of prediction of the temperature of the hottest spot in power transformer by solving the heat transfer partial differential equation (PDE) numerically. We have developed a numerical method based on the heat transfer theory using the finite element method and we only need to solve the heat conduction equation. The purely numerical approach for evaluating the hottest spot temperature and its location presented in this paper appears to provide results that correspond reasonably well with results of calculations and actual tests and on site measurements [9,10]. The authors are currently working on 3D model of transformer for estimation of temperature at different point of transformer knowing the load conditions. References [1] Pierce LW. Predicting liquid filled transformer-loading capability. IEEE Trans Ind Appl 1994;30(1):170–8. [2] IEC Standard, IEC60076-7. Loading guide for oil-immersed transformers; 2006.

[3] IEEE Standards, C57.91-1995. IEEE guide for loading mineral-oil-immersed transformer, vols. 18–19; 1996. p. 46–53. [4] IEEE Standards, 1538-2000. IEEE guide for determination of maximum winding temperature rise in liquid-filled transformers; 2000. [5] Radakovic Z. Numerical determination of characteristic temperatures in directly loaded power oil transformer. Eur Trans Electr Power 2003;13(1): 47–54. [6] Radakovic Z, Kalic DJ. Results of a novel algorithm for the calculation of the characteristic temperatures in power oil transformers. Electr Eng 1997;80(3):205–14. [7] Swift G, Molinski T, Lehn W, Bray R. A fundamental approach to transformer thermal modeling – Part I: theory and equivalent circuit. IEEE Trans Power Deliv 2001;16(2):171–5. [8] Tang WH, Wu OH, Richardson ZJ. Equivalent heat circuit based power transformer thermal model. IEEE Proc Electr Power Appl 2002;149. 2. [9] Pierce LW. An investigation of the thermal performance of an oil filled transformer winding. IEEE Trans Power Deliv 1992;7(3):1347–58. [10] Pradhan MK, Ramu TS. Prediction of hottest spot temperature (HST) in power and station transformers. IEEE Trans Power Deliv 2003;18(4): 1275–83. [11] Ryder SA. A simple method for calculating winding temperature gradient in power transformers. IEEE Trans Power Deliv 2002;17(4):977–82. [12] Swift G, Molinski TS, Bray R, Menzies R. A fundamental approach to transformer thermal modeling – Part II: field verification. IEEE Trans Power Deliv 2001;16(2):176–80. [13] Radakovic Z, Feser K. A new method for the calculation of the hot-spot temperature in power transformers with ONAN cooling. IEEE Trans Power Deliv 2003;18(4):1–9. [14] Reddy AS, Vijaykumar M. Hottest spot and life evaluation of power transformer design using finite element method. J Theor Appl Inform Technol 2008;4(3):238–43. [15] Faiz J, Sharifian MBB, Fakhri A. Two-dimensional finite element thermal modeling of an oil-immersed transformer. Eur Trans Electr Power 2007;18(6): 577–94. [16] Firouzifar S, Mahmoudi J. Assessment of power transformer cooler with FEMLAB. In: Proceedings of the 48th Scandinavian conference on simulation and modeling (SIMS 2007), 30–31 October. Göteborg (Särö), Sweden; 2007. [17] Mousavi Takami K, Mahmoudi J. Numerical modeling of heat generation and distribution in the core and winding of power transformers. Int J Emer Electr Power Syst 2008;9(2) [Article 7]. [18] Incropera FP, DeWitt DP. Fundamentals of heat and mass transfer. 4th ed. New York: J. Wiley and Sons; 1996.