- Email: [email protected]
Computation of temperature distributions in transformer covers due to high crossing currents in bushing regions
Electrical Power and Energy Systems 113 (2019) 699–712
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Computation of temperature distributions in transformer covers due to high crossing currents in bushing regions D. Cahue-Diaza, S. Maximova,
T
⁎,1
, R. Escarela-Perezb, J.C. Olivares-Galvanb, J. Alvarez-Ramirezc
a
Instituto Tecnologico de Morelia, PGIIE, Av. Tecnologico No. 1500, Lomas de Santiaguito, 58120 Morelia, Mich., Mexico Departamento de Energia, Universidad Autonoma Metropolitana-Azcapotzalco, 02200 Mexico City, Mexico c Departamento de Ingenieria de Procesos e Hidraulica, Universidad Autonoma Metropolitana-Iztapalapa, Mexico City 09340, Mexico b
ARTICLE INFO
ABSTRACT
Keywords: Power transformer Temperature distribution Analytical methods Heat equation Finite element method
Analytical formulae are widely applied in transformer design due to their easy implementation and low demand of computational resources. One of the great challenges in the field of transformer design is the acquisition of closed-form formulae for calculating temperature distributions in transformer tanks, taking account of different types of surrounding media. The presence of these media results in a variety of boundary conditions when solving the heat equation. In this research, a new rigorous analytical formula for calculating temperature profiles in bushing regions of transformer tank walls is obtained by solving the heat equation, with boundary conditions that properly model the heat interchange between different surrounding media. The heat source is modeled using a formula taken from a previously published investigation performed by the authors. Several study cases are performed, where analytical results are compared with finite element simulations for different boundary conditions. Analytical and numerical results show a good match, demonstrating high precision of our developed analytical method. Thus, this method provides a powerful tool for calculating temperature distributions in transformer tanks, which can be employed for optimization purposes during transformer design stages.
1. Introduction Power losses, that are originated in windings, core, tank covers, etc., result in heating of transformers. Overheating and presence of hot-spots can damage vital parts of a transformer, affecting its optimal operating regime and resulting in possible failures. Therefore, in order to improve transformer design and prevent undesirable consequences of such critical situations, it is vital to have means to compute temperature distributions and finding such hot-spots [1–5]. Several techniques for computing temperature distributions in power transformers have been developed throughout the history of transformer manufacture. The commonly used Finite Element Method (FEM) is a powerful and reliable tool for computing temperature profiles in active parts of transformers [5–9]. In particular, temperature distributions in tank walls [5,10] and oil temperature in power transformers [11–15] can be estimated by this method. In [6], FEM was employed to calculate the temperature distribution in a dry-type transformer. Joule losses in windings and eddy-current losses in magnetic cores were experimentally assessed. The results were employed to find hot-spots and air/oil flows inside the tank. The
authors emphasized the high accuracy level of the proposed experimental approach. In [7], FEM permitted the study of heat transfer and oil flow distribution, inside a 3-phase power transformer, by solving the continuity, momentum and energy balance equations in steady state. Different geometries and flow rates were analyzed in order to optimize the transformer cooling. In [8], fluid flow and convection heat transfer, in a thermo-convective circuit with several cooling channels and different inlet fluid velocities, were studied in a 3-phase power transformer immersed in mineral oil. The optimal velocity, that keeps temperature within stipulated limits, was obtained. In [9], the development of an advanced 3D finite element model, for the coupled solution of the heat transfer and fluid flow equations, that govern the transformer thermal performance, was presented. One of the main advantages of the proposed method is that there is no need to predefine convection coefficients at the interfaces between the active tank walls and circulating oil, since they are computed while the coupled heat-magnetic problem is being solved. In [12,14,15], numerical and experimental studies on a radiator-fan assembly of a power transformer were performed for different cooling configurations, to determine their cooling capacity and optimize their design.
Corresponding author. E-mail addresses: [email protected], [email protected] (S. Maximov). 1 On sabbatical leave: Universidad Autonoma Metropolitana-Azcapotzalco, 02200 Mexico City, Mexico. ⁎
https://doi.org/10.1016/j.ijepes.2019.06.019 Received 26 January 2019; Received in revised form 4 April 2019; Accepted 6 June 2019 Available online 14 June 2019 0142-0615/ © 2019 Elsevier Ltd. All rights reserved.
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Nomenclature
q q cr t kt T Ta W r r
convection coefficient global convection coefficient symmetric component of the temperature asymmetric component of the temperature Kronecker delta internal disk radius external disk radius metallic plate thickness steel electric conductivity steel relative permeability frequency
c
hc T+ T
heat flux by conduction heat flux by convection time thermal conductivity temperature ambient temperature volumetric heat source radiation emissivity Stefan-Boltzmann constant
nm
a b h
µr f
However, FEM requires considerable time-consuming computational resources, especially when dealing with 3D geometries and very pronounced skin-effects, where a high density meshing is normally required. Although numerical methods provide valuable information about thermal characteristics of transformers under fixed conditions (geometry, boundary conditions, etc.), they cannot provide information of how the temperature distribution depends explicitly on each parameter of the transformer. On the other hand, analytical formulae explicitly contain all the parameters of the transformer and surrounding media. This fact makes analytical methods powerful tools for design optimization of transformers using its parameters (geometry, materials involved in the transformer construction, etc.). Another aspect, where analytical approaches are useful, is involved with the determination of heat transfer coefficients. It is well-known that analytical or experimental determination of these coefficients is a difficult task [16,4]. Once an analytical solution to the heat equation is obtained, it can be compared with measurements. Analytical models can be fitted to experimental points by variation of heat transfer coefficients in the analytical solution, using least squares [16] or multiobjective (MO) optimization procedures [4]. Heat transfer coefficient values, thus obtained, can be consequently employed for transformer design purposes. Recently, an analytical formula to calculate temperature distribution in the bushing region of transformers has been obtained [16]. The heat source, employed to calculate the temperature profile, can be considered numerically or analytically. The proposed formula agrees well with FEM calculations. Nevertheless, it was assumed in [16] that the tank wall was surrounded by an homogeneous medium (e.g., air) with the same ambient temperature on each side and, therefore, the same heat transfer coefficient throughout the boundary between the metallic plate and the medium. Hence, the assumptions that were adopted in that work are far from real. In this paper, a new mathematical model for calculating temperature profiles in bushing regions of power transformers, that takes account of different media on each side of the tank wall, along with their heat transfer mechanisms (through boundaries), is derived. The influence of media on different sides of the wall is achieved by means of ambient temperature and heat transfer coefficients through the interface. The heat transfer coefficients depend on three mechanisms: conduction, convection and radiation. Several study cases are carried out in order to show that a good match between the new formula and FEM results is achieved for different values of the ambient temperature inside the tank. Here, the implementation of the analytical solution in the determination of the heat transfer coefficients through MO optimization procedure is discussed.
tank cover producing a nonuniform temperature distribution in the bushing region as a result of three processes: Ohm heating, thermal conduction and cooling through the surface. The tank geometry in the bushing region is presented in Fig. 1. Fig. 1B shows the cross-section of the cover of the transformer that passes along the conductor. The schematic view of this cross-section is presented in Fig. 1C. Fig. 2 shows the boundaries through which the transformer tank cover interchanges heat with different surrounding media. The upper surface 3 of the plate is in contact with air, so that the two heat transfer mechanisms through this boundary, namely convection and radiation, are to be considered. 4 is the inner surface of the transformer tank, through which the tank wall interacts with the internal medium of the transformer. The tank can be filled with air or oil, so that the thermal interaction of the tank cover with the fluid inside the transformer is given by convection heat transfer. The space inside the hole, which is a dielectric, comes into contact with the tank wall through surface 1, whereas the medium on the side of boundary 2 is air. The heat equation that governs the heat transfer phenomenon in steady state (see [16]) is given by
1 T r r r r
+
2T
z2
=
1 W r, z , kt
(1)
where the time-averaged heat source W (r , z ) can be calculated numerically [17] or analytically [18]. Due to the symmetry of the problem, cylindrical coordinates are employed. Boundary conditions on each surface can be expressed as follows:
nk ·q cr
k
= nk · q
k,
where q is the heat flux density in the tank plate, which is expressed by
2. The model A transformer bushing is an insulator that allows an electrical conductor to pass safely through the tank wall. Alternating current in the conductor produces an electromagnetic field around it, which, in turn, generates eddy currents in the tank wall. These currents heat the
Fig. 1. Bushing region of a transformer. 700
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
kt
Fourier’s law as:
kt T ,
c (T )(T
Ta) +
r
r
(T 4
Ta4 )} n ,
q cr =
(T ) +
r
T 3 mTam
r
T
T± r , z =
Ta n.
The following heat transfer coefficient is therefore obtained:
1 r r
3
(T ) +
r
T3
r
mT m a
m =0
and can be considered as a constant parameter for a given range of temperatures. Then, we find that
q cr = hc (T
Ta ) n.
kt
kt
h c1 ( T
Ta1 )
= 0, r=a
T + hc 2 (T r
Ta2 )
T + h c 3 (T z
Ta3 )
= 0, r=b
= 0, z = h /2
T (r , z ) ± T (r , 2
(r ) + T+ r
z)
1 W kt
=
T+ r
hc1 (T+
Ta1 )
kt
T+ r
+ hc 2 (T+
Ta2 )
+ h+ (T+
Ta +)
1 r r
(r ) + T r
2T
T r
hc1 T
T kt r
+ hc 2 T
kt
kt
(3)
2T + z2
kt
T kt z+
Experimental determination of the average parameter hc is discussed in Appendix B. Now, the four boundary conditions take the following form:
T kt r
(6)
.
(7)
Since the heat source W (r , z ) is an even function of z (see [18]), Eq. (1), along with boundary conditions (3)–(6), can be divided into the following two problems:
m=0
h c (T ) =
= 0, z = h /2
where T+ (r , z ) and T (r , z ) are the symmetric and antisymmetric components (with respect to the coordinate z) of the temperature distribution, which are defined in the following way:
(2)
where r is the Stephan-Boltzmann constant and n is a unitary vector normal to the surface. In [16], a linearization procedure for the convection heat flux density has been introduced. According to this method, Eq. (2) should firstly be cast in the form: 3
Ta 4 )
T (r , z ) = T+ (r , z ) + T (r , z ),
and q cr consists of two convection and radiation heat exchange terms:
q cr = {
h c 4 (T
where the existence of four different heat transfer coefficients on each surface and four different ambient temperatures near each side of the tank wall has been assumed. In order to solve Eq. (1) with the asymmetric boundary conditions (3)–(6), it is appropriate to divide the temperature field T (r , z ) into two parts as follows:
Fig. 2. Tank geometry.
q=
T z
T z
z2
+ h+T
r=a r=b
r, z , = 0, = 0,
z = h /2
=
h T |z = h /2 ,
(8)
= 0,
r=a r=b
= 0, = 0,
z = h /2
=
h (T+
Ta )z = h /2 ,
(9)
where (4)
(5)
h± =
h c3 ± h c 4 , 2
(10)
Ta ± =
hc 3 Ta3 ± hc 4 Ta 4 . h c3 ± h c 4
(11) Fig. 3. Schematic behavior of the left-hand part of (29) and solutions for qn .
701
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
where T±[N ] (r , z ) is the N th iteration for the temperature T± (r , z ) . Here, the zero iteration for the antisymmetric part of temperature distribution is taken as T [0] (r , z ) 0 . The scheme (13) and (14) can be truncated at a given iteration that provides sufficient precision of the approximate solution. The hypothesis of this work is that, with a good precision, the iteration scheme given by (13) and (14) can be considered up to the first iteration. Further, this hypothesis is verified by comparing analytical results with numerical computations.
Problems (8) and (9) are coupled through boundary conditions, so they cannot be solved separately. In the next section, the scheme of solution of these two coupled problems is developed. 3. Proposed iteration scheme Problems (8) and (9) can be decoupled if the right-hand part of the last equation of problem (8) is neglected. This term can be estimated using the mean value theorem, according to which
z
h/2, h/2 :
T r, z =
T (r , z )
T (r ,
z)
2
T = z
3.1. First iteration for the symmetric part
z.
The problem for the first iteration of the symmetric part of the temperature distribution T+ (r , z ) can be obtained from (13) by substituting N = 1 and T [0] (r , z ) 0 . The problem thus obtained displays similarity to the problem solved in [16]. However, in this case different heat transfer coefficients hc1, hc 2 and ambient temperatures Ta1 and Ta2 near the boundaries 1 and 2 are now accounted for. Following the mathematical steps given in [16], we obtain:
z=z
The following estimation can therefore be demonstrated:
h T r, z
h
T h h h= z kt
max
h /2 z h /2
max
h /2 z h /2
qz ,
(12)
where max h /2 z h /2 qz is the maximum heat flux density in the z-axis direction. Hence, this term is small if the wall thickness h is also small enough and the difference of ambient temperatures above and below the tank wall is not large. The last statement means that the heat flux density in the z-direction is small (Fourier’s law). Nevertheless, when overheating of the medium below the cover plate takes place, the right-hand part of the last equation of problem (8) becomes significant and, consequently, it cannot be neglected. In this case, an iteration scheme can be applied. In the first iteration, the righthand part of the last equation of (8) is neglected. As a result, problem (8) is solved independently of (9). In the next step, the solution of this problem is substituted into (9). In the second iteration, the solution to (9) is substituted into (8), where the right-hand part of boundary conditions is now taken into account. These consecutive steps can be represented by the following iteration scheme: 1 r r
kt kt kt
1 r r
r
[N ] T+
r
[N ] T+
2T [N ] + z2
hc1 (T+[N ]
r
1 W kt
=
= Ta + + Wn ( ) d
r
+ hc 2 (T+[N ]
z
+ h+ (T+[N ]
= 0, =
+
2T [N ]
z2
kt
T [N ] r
hc1 T [N ]
kt
T [N ] r
+ hc 2 T [N ]
T [N ] kt z
+ h+T
h T [N
1]
r , h/2 ,
(13)
= 0,
+
1 kt
b a
Gn r ,
Dn(12) I0 ( +
(11) n ) + Dn K 0 ( n (11) (22) Dn Dn Dn(12) Dn(21) (12) Dn hc 2 a2 h sin n , 2 Dn(12) Dn(21)
)
Wn ( ) d
n
2 Zn
2
(16)
Dn(11) kt
b
Dn(22) I0 (
a a1
+
(21) n ) + Dn K 0 ( n (11) (22) Dn Dn Dn(12) Dn(21) Dn(11) hc 2 a2 h sin n , 2 Dn(12) Dn(21)
Dn(11) = kt
n I1 ( n a)
Dn(12) Dn(21) Dn(22)
= kt
n K1 ( na )
= kt
n I1 ( n b)
= kt
n K1 ( nb)
)
Wn ( ) d
n
2 Zn
2
(17)
hc1 I0 ( na), + hc1 K 0 ( na),
+ hc 2 I0 ( nb), (18)
hc 2 K 0 ( nb).
Here,
= 0,
Gn (r , ) = I0 ( nr ) K 0 (
r=b
=
b
Dn(11) Dn(22)
r=a
z = h /2
nr
(15)
a
Dn(21) hc1
= 0,
[N ]
+ Bn K 0
Bn =
Ta +)
nr
,
Dn(22) hc1 a1 Dn(11) Dn(22)
z = h /2
T [N ] r
Dn(22) kt
=
r=b
[N ] T+
An I0
An
= 0,
Ta2 )
nz
where
r, z ,
Ta1 )
cos n=0
r=a
[N ] T+
r
+
T+[1] r , z
h
(T+[N ] (r ,
h/2)
Ta ),
n
) (
r ) + I0 (
n
) K 0 ( nr )[1
(
r )] (19)
(14)
is the Green’s function, where I0 ( nr ) and K 0 ( nr ) are modified Bessel
Fig. 4. Mesh used for FEM analysis. 702
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Finally, {
Table 1 Parameters for simulations. Parameters
Case I
Case II
Case III
µr kt f h c1 h c2 h c3 h c4 Ta1 Ta2 Ta3 Ta4
500 16 W/mK 50 s−1 25 W/m2K 10 W/W/m2K 10 W/W/m2K 35 W/W/m2K 80 °C 18 °C 18 °C 80 °C
500 16 W/mK 50 s−1 25 W/W/m2K 10 W/W/m2K 10 W/W/m2K 35 W/W/m2K 100 °C 18 °C 18 °C 100 °C
500 52 W/mK 50 s−1 25 W/W/m2K 10 W/W/m2K 10 W/W/m2K 35 W/W/m2K 100 °C 18 °C 18 °C 100 °C
n
{0, }} are solutions to the following equation:
n
h h = +. 2 kt
tan
(21)
Eq. (21) can be solved either numerically or analytically (see [16]). Here, the asymptotic analytical solution to Eq. (21), obtained in [16], is employed. 3.2. First iteration for the antisymmetric component Solution (15) should be substituted into the right-hand side of the last equation of (9) with N = 1. The solution to (9) can be sought using separation of variables [16] as follows: (22)
T [1] (r , z ) = R (r ) Z (z ). As a result, the following equations can be obtained:
R +
1 R + q2R = 0, r
(23)
Z
q2 Z = 0
(24)
with boundary conditions
[kt R hc1 R]r = a = 0, [kt R + hc 2 R]r = b = 0.
(25)
General solutions to Eqs. (23) and (24) can be obtained in the form:
Wn (r ) =
2 h + sin( nh)/
h /2 n
h /2
nz
dz.
Z = sinh(qz ),
(27)
(kt qJ1 (qa) + hc1 J0 (qa)) C1 + (kt qY1 (qa) + hc1 Y0 (qa)) C2 = 0,
r ) is Heaviside unit-step function.
W r , z cos
(26)
where J0 (qr ) and Y0 (qr ) are Bessel functions of zero order and q is an unknown constant. Substitution of the general solution (26) into the boundary conditions (25) yields the following linear system of algebraic equations:
Fig. 5. 2D temperature distribution. FEM-simulations. Case I.
functions of zero order, while (
R = C1 J0 (qr ) + C2 Y0 (qr ),
(kt qJ1 (qb)
hc 2 J0 (qb)) C1 + (kt qY1 (qb)
hc 2 Y0 (qb)) C2 = 0.
(28)
The existence of a nontrivial solution to system (28) is provided if the determinant of the coefficient matrix of system (28) is equal to zero:
(20)
Fig. 6. Temperature profiles on the upper surface: theoretical and FEM-calculated curves. Case I. 703
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Fig. 7. Temperature profiles on the lower surface: theoretical and FEM-calculated curves. Case I. n (r )
Det = (kt qJ1 (qa) + hc1 J0 (qa))(kt qY1 (qb) (kt qJ1 (qb)
hc 2 Y0 (qb))
hc 2 J0 (qb))(kt qY1 (qa) + hc1 Y0 (qa)) = 0,
= (kt qn Y1 (qn a) + hc1 Y0 (qn a)) J0 (qn r ) (29)
} are orthogonal, as shown in Appendix A. The functions { n (r ) n The norm of these functions is also calculated there. Therefore, instead of this system, the following set of orthogonal and unitary functions can be considered:
which is the base to find the values of constant q. Eq. (29) has an in} , which can be numerically obtained. finite number of roots {qn n The solutions for q are shown in the Fig. 3. After substituting qn into (28) and taking account of (29), the following general solution to system (28) is obtained:
C1, n = Cn (kt qn Y1 (qn a) + hc1 Y0 (qn a)) C2, n = Cn (kt qn J1 (qn a) + hc1 J0 (qn a)).
n (r )
n (r )
=
n (r )
.
(32)
n
As a result, the general solution for the asymmetric part of the temperature distribution takes the form:
(30)
Thus, the solutions to Eqs. (23) and (24) take the form:
Rn (r ) = Cn
(kt qn J1 (qn a) + hc1 J0 (qn a)) Y0 (31)
(qn r ).
T [1] r , z =
and Zn (z ) = sinh(qn z ),
Cn
n (r )sinh
qn z .
n=1
(33)
Substitution of this solution into the last boundary condition of problem
where
Fig. 8. Temperature profiles on the upper and lower surfaces. Analytical solution. Case I. 704
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
T [1] r , z h sinh(qn z )
= T+[1] r , z + n= 1
×
(9) yields:
Cn
n (r )(qn kt cosh(qn h /2)
=
h (T+[1] (r , h/2)
+ h+sinh(qn h/2))
h
b a
(Ta
T+[1] ( , h /2))
n(
) d
qn kt cosh(qn h/2) + h+sinh(qn h /2)
.
T+[1] ( , h/2))
n(
) d .
(36)
Three different study cases are considered in this section to illustrate the accuracy of our analytical approach. Ansys Fluent software has been employed to carry out FEM simulations considering an axisymmetric geometry (see Fig. 2). A cylindrical metallic plate with internal and external radii a = 0.0125 m and b = 0.265 m, respectively, and a thickness h = 0.01 m, has been considered. These geometry parameters were used for all the simulations. The cross-section of the tank wall was covered using a mesh with a refinement near the plate surface, as shown in Fig. 4, so as to accurately compute eddy current losses (the power density of the eddy current losses is higher near the tank wall surface [18,19]). The electromagnetic model reported in [18] was employed for calculating eddycurrent losses while our new analytical method is used to compute the temperature distribution. Since the model [18] faithfully reproduces the losses, there was no need to consider different currents in the
Using the orthogonality property of the system of functions { n (r ) n } , the following expression for the unknown constants Cn can be obtained from Eq. (34):
Cn =
(Ta
4. Study cases (34)
Ta ).
a
Eq. (36), together with (15), represents a new model for calculating temperature distributions in bushing areas of transformer tanks. The implementation of this model requires: the heat source function W (r , z ) (calculated numerically or analytically), the geometry parameters, properties of the tank wall material and surrounding fluid and the heat transfer coefficients on the tank wall surface. The algorithm for calculating the temperature distribution consists of the following steps: (i) calculation of average ambient temperatures (11) and heat transfer coefficient (10) on surfaces 3 and 4 ; (ii) estimation of the parameters n by solving Eq. (21) numerically or analytically (see [16]); (iii) computation of Fourier coefficients (20) of the heat source; (iv) calculation of the parameters (18) and constants (16) and (17); and (v) finally, implementation of formula (36). Finally, the appropriate number of terms that should be taken from series (36), to get a good precision, is estimated in Appendix C.
Fig. 9. 2D temperature distribution. FEM-simulations. Case II.
n=1
n (r )
b
qn kt cosh(qn h /2) + h+sinh(qn h/2)
(35)
Substitution of (35) into (33) and summing the symmetric (15) and antisymmetric (33) components of the solution, the following temperature distribution in the first order iteration is obtained:
Fig. 10. Temperature profiles on the upper surface: theoretical and FEM-calculated curves. Case II. 705
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Fig. 11. Temperature profiles on the lower surface: theoretical and FEM-calculated curves. Case II.
conductor to represent different levels of power loss. Therefore, all the simulations were performed using the same value I = 1000 A of the current in the conductor, i.e., the same volumetric heat source was considered in each case. So, the simulations of our new work were focused on modeling different boundary conditions round the tank wall, namely, different ambient temperatures. The operating conditions and parameter values of the system are given in Table 1. The loss density employed for the analytical computations of the temperature distributions was calculated using the electromagnetic model published by the authors in [18]. Numerical computations were performed using FEM both for calculating power losses and for computing the temperature distribution. The results of the simulations are presented in Figs. 5–16. Case I corresponds to the temperature of the medium below the transformer cover Ta4 = 80 °C and thermal conductivity of the tank wall kt = 16 W/mK . Fig. 5 shows the axisymmetric temperature distribution obtained on the upper surface of the cover, using FEM-simulations. Figs. 6 and 7 compare the results of FEM-simulations with our
analytical results, obtained with formula (36), on the upper and lower surfaces of the tank wall, respectively. It can be observed in Figs. 6 and 7 that both numerical and analytical results closely match. In Fig. 8, the difference between analytically calculated temperatures on the upper and lower surfaces, due to difference between ambient temperatures above and below the tank wall, can be observed. In the second case, the temperatures below the covering plate and in the hole were set to Ta4 = 100 °C and Ta1 = 100 °C , in order to model oil overheating in the transformer and providing a higher temperature gradient in the tank wall thickness. The numerical and analytical results are presented in Figs. 9–12. It can be observed in these figures that the increment of 20 °C in the temperature near the lower surface leads to increased temperature values at each point of the tank wall. However, it did not noticeably affect the temperature gradient in the tank wall thickness. The third case (Figs. 13–16) showed that the increment of the thermal conductivity of the tank wall from 16 W/mK to 52 W/mK considerably decreases the temperature gradient in both the axial and
Fig. 12. Temperature profiles on the upper and lower surfaces. Analytical solution. Case II. 706
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
exchange through surface 2 is not an important factor, since the area of this surface is small with respect to 3 and 4 (the plate thickness h is small with respect to the tank dimensions). Moreover, the heat-exchange velocity near 2 is much smaller than in the bushing region, since the temperature near the outer border is smaller. Therefore, the surface 2 geometry can be disregarded. However, the total area of boundaries 3 and 4 is an important factor in the heat transfer problem. In this section, it is shown that our analytical development can be generalized for actual geometries such as that presented in Fig. 1, where the bushing is displaced from the center of the plate. It is assumed that the radius of the area, where the main part of losses is concentrated (e.g., 90% of total losses), is much smaller than the outer dimensions of the cover. To do so, we consider formulae (15)–(19). In virtue of the asymptotic behavior of the modified Bessel functions (Kn ( nr ) rapidly decreases whereas In ( nr ) increases as the distance r grows [20]), it follows from (18) that 1/ Dn(21) and Dn(22) are small parameters. Formally, 0 . Application of this fact to Eqs. it means that Dn(21) and Dn(22) (16) and (17) shows that the parameter An is small enough to be neglected in Eq. (15) and the parameter Bn can be simplified as follows: Fig. 13. 2D temperature distribution. FEM-simulations. Case III.
Bn
radial directions. The three analyzed cases demonstrate that our analytical formula (36) is reliable for calculating temperature distributions and evaluating hot spots in bushing regions of transformer tank walls. The selected iteration scheme (13) and (14), truncated at the first order iteration resulted in a good precision when computing temperature distribution even in the case of an extremely high temperature below the tank wall.
1 Dn(11) kt Dn(12)
b a
K0
Wn ( ) d +
n
n
2 Zn
2
hc1 a1 sin Dn(12)
nh . 2
(37) As a result, solution (15) takes the form:
T+[1] r , z = Ta + +
cos
nz
{
n=0
5. Generalization for nonaxisymmetric geometries
+
Results (15) and (36) were obtained for a geometry with axial symmetry. However, bushings are mounted at points displaced from the center of the tank cover in many cases. Therefore, applicability of the obtained results for real cases, such as that shown in Fig. 1A, should be verified. Otherwise, solution (36) should be generalized for other geometries, the path that is undertaken next. It is well-known that the heat source is mainly concentrated in an area near the conductor, so that the main contribution to the tank wall heating is localized in the bushing region. It is also clear that the heat
n
2 Zn
2
h c1 a 1 sin Dn(12)
1 Dn(11) kt Dn(12) nh
2
K0
b a
K0
nr
Wn ( ) d
n
+
1 kt
b a
Gn r ,
Wn ( ) d }. (38)
The analysis of this formula shows that: (a) formula (38) does not depend on the ambient temperature T2 , on the right side of the covering plate; (b) the dependence of (38) on the outer shape of the transformer cover is only given by the upper limit of the integrals in the right-hand side of (38). Hence, the dependence of (38) on the covering plate shape is very weak. It can be proved by taking a derivative of (38) with
Fig. 14. Temperature profiles on the upper surface: theoretical and FEM-calculated curves. Case III. 707
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Fig. 15. Temperature profiles on the lower surface: theoretical and FEM-calculated curves. Case III.
respect to the outer radius b that
T+[1] = b
cos n=0
nz
{
1 Dn(11) K0 kt Dn(12)
where Eq. (40) represents a convex curve. Thus, solution (41) can be generalized as follows: nr
K0
nb
+
1 Gn r , b } Wn (b) b, kt
T+[1] r , z
(39) which is negligible, since the heat source W (r , z ) (and consequently the Fourier coefficients Wn (r ) ) decreases rapidly by increasing the distance r. Moreover, the values of the Bessel function K 0 ( nb) on the cover border b is also negligible. This and the fact that Eq. (38) does not depend on the temperature near the boundary 2 and the convection coefficient hc2 , proves that solution (38) depends weakly on the boundary 2 shape. In turn, this allows the generalization of solution (38) for the case of a transformer cover of an arbitrary geometry. Let us consider a transformer cover with a shape given by the following equation (see Fig. 17):
r = b ( ),
[0, 2 ],
= Ta + +
cos
nz
{
n=0
+ n b( ) a
2 Zn
2
h c1 a 1 sin Dn(12)
Gn r ,
Dn(11) 1 Dn(12) 2 kt nh 2
K0
2 0
nr
Wn ( ) d },
d
+
b( ) a
1 2 kt
K0 2 0
n
Wn ( ) d
d
(41)
where the integrals in the right-hand side are integrated over the boundaries 3 and 4 . In the case of b = const , Eq. (41) is reduced to (38).
(40)
Fig. 16. Temperature profiles on the upper and lower surfaces. Analytical solution. Case III. 708
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
result, we obtain:
T [1] r , , z =
T [N ] 1 1 r + 2 r r r r
2T [1] 2
+
2T [N ]
z2
T [N ] 1 r + r r r
2T [N ]
z2
= 0,
which proves that generalization (42) satisfies the heat equation even for the nonaxisymmetric case. In conclusion, it should be noticed that b = const means that Eq. (42) takes the form (36). Thus, a good precision solution, generalized for the case of an arbitrary geometry, takes the form:
T [1] (r , , z ) = T+[1] (r , , z ) + T [1] (r , , z ), where T+[1] (r , , z ) and T [1] (r , , z ) are represented by (41) and (42), respectively. 6. Conclusion Fig. 17. Asymmetric geometry of the transformer cover.
In this paper, a new analytical formula for calculating temperature distributions in bushing regions of transformer tanks has been obtained. The new formula takes account of the fact that the tank wall could be surrounded by different media with different thermal properties such as the heat transfer coefficient, ambient temperature and heat radiation emissivity. In particular, the plate can be in contact with cooling fluids on the inner side (e.g., oil) and air on the outer side. Unlike previous works, all these possibilities have been taken into account in our proposed formula, through lumped parameters of heat transfer coefficients and different ambient temperatures. As a result, the new analytical approach is capable to reproduce faithfully temperature gradients in the tank wall thickness, produced by the difference of ambient temperatures above and below the tank wall, which was impossible with previous analytical models (see, e.g., [16,4]). The analytical results have been compared with FEM simulations, which have also been performed considering heat source calculations based on electromagnetic FEM calculations. The comparison shows that our new formula for temperature distributions leads to a good match with numerical results. In turn, this demonstrates that the proposed analytical method is reliable for predicting temperature distributions in transformer tanks and estimating hot spots. The analytical formula provides high precision results with no need of expensive computational resources. Unlike numerical methods, the analytical formula explicitly depends on all the parameters of the system such as geometry and material parameters involved in the transformer construction. This makes our analytical method preferable for transformer optimal design with respect to its parameters and for determination of thermal parameters through MO optimization procedures.
A similar arbitrary geometry generalization for the antisymmetric part of (36) can be realized as follows:
T [1] r , , z
n=1
h sinh(qn ( ) z ) × qn ( ) kt cosh(qn ( ) h/2) + h+sinh(qn ( ) h/2)
1 2
0
=
2
d
b( ) a
(Ta
T+[1] ( , h/2))
,
n
n
r,
d .
(42)
Here, the dependence of the parameters qn on the angle can be obtained from (29), where parameter b is a function of (see (40)). Finally, functions n (r , ) , as well as Eq. (42), depend on the angle through parameters qn and b. Let us show that (42) satisfies the heat equation with good precision. To do so, the derivative of (42) with respect to should be estimated:
T [1]
T [1] qn T [1] b + = qn b
=
T [1] dqn T [1] n db + . qn db d n db
(43)
In order to estimate the derivative of dqn / db , a derivative of Eq. (29) with respect to b should be taken. As a result, after some bulky, but easy calculations, the following estimation can be obtained:
dqn db
qn
=
+
b
1 , b2
(44)
where is the Landau big O notation. In other words, parameters qn depend weakly on the cover shape. In a similar way, the estimation of the derivative of n (r ) with respect to b shows that it is negligible. Summarizing, we are led to the conclusion that derivative (43) is small
(1/b 2)
enough. Similarly, it can be shown that
2T [1] 2
Acknowledgements The authors are grateful for financial support provided by the following Mexican Institutions: CONACYT, SNI and PRODEP.
is also negligible. As a
Appendix A
} can be proved in a standard way. Let us consider two functions of this set: The orthogonality of the set of functions { n (r ) n These functions satisfy Eq. (23) for different values of q and boundary conditions (25): m
+
1 r
m
n
+
1 r
n
where qm
+ qm2
+ qn2
m
n
= 0,
and
n (r ) .
(45)
= 0,
qn . Multiplication of Eq. (45) by r
m (r )
(46) n (r )
and Eq. (46) by r
m (r )
709
and subtraction of the second equation from the first one, leads, after
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
simplifying and integrating with respect to r over the segment [a , b], to the following equation:
r(
n
m
b n )|a
m
+ (qm2
b
qn2 )
a
m
n rdr .
(47)
Substitution of boundary conditions (25) into (47) leads to the following orthogonality condition b
n rdr
m
a
=0
n. for m In order to find the norm of functions (q, r ) = (kt qn Y1 (qn a) + hc1 Y0 (qn a)) J0 (qr ) instead of m , where q can be obtained: n
b
2 n rdr
a
=
r(
lim
q
n q2
qn
qn2
let us substitute into Eq. (47) the function (48)
(kt qn J1 (qn a) + hc1 J0 (qn a)) Y0 (qr )
qn . It is evident that
2
=
n (r ) ,
b n )|a
=
(q, r ) satisfies Eq. (23) and
1 lim r 2qn q qn
(qn , r ) =
n (r ) .
Then, the following expression for the squared norm
b
q
n
q
,
n
(49)
a
where L’Hôpital’s rule has been used to calculate the limit. From (48), it can be directly verified that
q
=
qr
and
=
q
r q
.
(50)
Substitution of the properties (50) into (49) yields: 2
=
n
r2 (q 2 2qn2 n
2 n
+
2 b n )|a .
This result, with boundary conditions (25) taken into account, finally takes the form: 2
=
n
h2 1 b2 qn2 + c22 2 2qn kt
2 n (b )
a2 qn2 +
hc21 kt2
2 n (a )
.
(51)
Appendix B The difficulty of the exact theoretical determination of thermal parameters such as the heat transfer coefficient hc , radiation emissivity, etc., is well known (see [4,5,16]). Moreover, convection coefficients are, in general, temperature-dependent functions. Also, a direct determination by experiment is often complicated [21]. However, temperature assessment of transformers is sensitive to the precision of these parameters. Since our new formula (36) showed a good match with FE results, it can be effectively used for calculating thermal parameters. Recently, multi-objective optimization was efficiently employed for thermal parameter determination [4] through comparison of FEM simulations with measurements. The precision of our new formula permits its implementation within MO optimization, where thermal parameters hc1, hc 2, hc3 and hc4 are employed as independent variables. Let us denote Tk as the temperature values measured at different points on the surface of the metallic disc, with coordinates rk , where 4 is the following vector: k = 1, 2, …, N , by N sensors, whereas analytical results at the same points, obtained from (36), are: T (rk, H) , where H
H = (hc1, hc 2 , hc3, hc 4 ). Thus, the first objective function to be minimized is established as the N
f1 (H) =
k=1
2 -space
error:
Tk ) 2 .
(T (rk, H)
(52)
The minimization of only one parameter (52) with respect to the vector H (the so-called single-objective optimization) is the widely-known leastsquare method. Another objective function that was proposed in [4] is the sensitive function that is determined as follows:
f2 (H) =
max1
i 4 f1i
f1
f1
,
(53)
where f1i is the first objective function (52) corresponding to ± 4% relative perturbations of the parameters hci , where i = 1, 4 . Thus, the 2D vector f = (f1 , f2 ) should be optimized with respect to the 4D vector H . The advantage of the use of the analytical solution instead of FEM simulations in determining thermal parameters is that the analytical calculations are precise and do not require high-consuming computational resources. Appendix C The right-hand side of (36) is a generalized Fourier series with orthogonal and unitary system of functions { Parseval’s identity, namely, by the convergence of the series
n (r )} .
Its convergence is provided by
2
cn
<
,
(54)
n= 1
where 710
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
cn =
h sinh(qn z ) qn kt cosh(qn h/2) + h+sinh(qn h/2)
(55)
are coefficients of the series in (36). The number of terms that should be taken into account in (36) depends on the required precision, and it can be estimated based on series (54). Firstly, the minimum number of terms N is so that qN is big enough. This fact allows the following estimation for coefficients (55):
h sinh(qn h /2)
cn
qn kt cosh(qn h/2) + h+sinh(qn h/2)
=
h qn kt coth(qn h/2) + h+
h . qn kt
(56)
Therefore, series (54) is majorized by another one as follows: 2
cn n= 1
n= 1
h2 qn2 kt2
(57)
Secondly, parameters qn can be asymptotically estimated from (29) for the high values of n. Substitution of the asymptotic formulae for Bessel functions [20] into Eq. (29), after neglecting small terms, yields:
sin(qn (b
a))
0.
The solutions to this equation are:
n
qn
b
a
,
where n 2
cn n= 1
(58) . Thus, the convergence of series (57) can be estimated in the following way:
h2 kt2
n=1
1 qn2
h 2 (b kt2
a) 2 2
n=1
1 , n2
where
n= 1
2 1 = . 2 n 6
(59)
Let us truncate series (59) at the
rN = n=N +1
1 = n2
N th
term. Then, the residual term rN of the truncated series is:
N+1 ,
(60)
where
N+1 =
d ln (x ) dx
x=N+1
is the digamma function [22]. Therefore, in order to provide a good precision, we can make the residual term do not exceed a small part of the 1. Therefore the minimum number N can be found from the following complete series value. Namely, let be a small parameter such that 0 < equation:
N+1 For example, if N 11.
2
6
.
= 0.01 (i.e., the relative error does not exceed 1% ), the approximate value of N is 60. In the case of
= 0.05 (5% of the relative error)
[8] Chereches NC, Chereches M, Miron L, Hudisteanu S. Numerical study of cooling solutions inside a power transformer. Energy Procedia 2017;112:314–21. [9] Tsili M, Amoiralis EI, Kladas AG, Souflaris AT. 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. [10] Sitar R, Šulc I, Janić Ž. Prediction of local temperature rise in power transformer tank by fem. Procedia Eng 2017;202:231–9. [11] Najar S, Tissier JF, Cauet S, Etien E. Improving thermal model for oil temperature estimation in power distribution transformers. Appl Therm Eng 2017;119:73–8. [12] Anishek S, Sony R, Kumar JJ, Kamath PM. Performance analysis and optimisation of an oil natural air natural power transformer radiator. Procedia Technol 2016;24:428–35. [13] Alegi GL, Black WZ. Real-time thermal model for an oil-immersed, forced-air cooled transformer. IEEE Trans Power Deliv 1990;5(2):991–9. [14] Paramane SB, Van Der Veken W, Sharma A. A coupled internal-external flow and conjugate heat transfer simulations and experiments on radiators of a transformer. Appl Therm Eng 2016;103:961–70. [15] Garelli L, Ríos Rodriguez G, Storti M, Granata D, Amadei M, Rossetti M. Reduced model for the thermo-fluid dynamic analysis of a power transformer radiator working in ONAF mode. Appl Therm Eng 2017;124:855–64. [16] Maximov S, Escarela-Perez R, Olivares-Galvan JC, Guzman J, Campero-Littlewood
References [1] Allahbakhshi M, Akbari M. Heat analysis of the power transformer bushings using the finite element method. Appl Therm Eng 2016;100:714–20. [2] Cinar MA, Alboyaci B, Sengul M. Comparison of power loss and magnetic flux distribution in octagonal wound transformer core configurations. J Electr Eng Technol 2014:1290–5. [3] Díaz-Chacón JM, Hernandez C, Arjona MA. Finite element and neural network approach for positioning a magnetic shunt on the tank wall of a transformer. IET Electr Power Appl 2016;10:827–33. [4] Penabad Durán P, Di Barba P, Lopez-Fernandez X, Turowski J. Electromagnetic and thermal parameter identification method for best prediction of temperature distribution on transformer tank covers. Int J Comput Math Electr Electron Eng 2015;34:485–95. [5] Penabad-Duran P, Lopez-Fernandez XM, Turowski J. 3D non-linear magnetothermal behavior on transformer covers. Electr Power Syst Res 2015;121:333–40. [6] Ortiz C, Skorek AW, Lavoie M, Bénard P. Parallel CFD analysis of conjugate heat transfer in a dry-type transformer. IEEE Trans Ind Appl 2009;45:1530–4. [7] Wakil NE, Chereches NC, Padet J. Numerical study of heat transfer and fluid flow in a power transformer. Int J Therm Sci 2006;45(6):615–26.
711
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al. E. New analytical formulas for temperature assessment on transformer tanks. IEEE Trans Power Deliv 2016;31(3):1122–31. [17] Olivares-Galvan JC, Perez RE, Kulkarni SV, Leon FD, Venegas-Vega MA. 2d finiteelement determination of tank wall losses in pad-mounted transformers. Electr Power Syst Res 2004;71:179–85. [18] Maximov S, Olivares-Galvan JC, Escarela-Perez R, Magdaleno-Adame S, CamperoLittlewood E. New analytical formulae for electromagnetic field and eddy current losses in bushing regions of transformers. IEEE Trans Magn 2015;51. 6 300 710–1–10.
[19] Maximov S, Olivares-Galvan JC, Magdaleno-Adame S, Escarela-Perez R, CamperoLittlewood E. Calculation of nonlinear electromagnetic fields in the steel wall vicinity of transformer bushings. IEEE Trans Magn 2015;51. 6 300 811–1–11. [20] Olver FW. Asymptotics and special functions. 3rd ed. Academic Press; 1974. [21] Kitak P, Glotic A, Ticar I. Real-time thermal model for an oil-immersed, forced-air cooled transformer. IEEE Trans Magn 2015;50(2):933–6. [22] Abramowitz M, Stegun IA. Handbook of mathematical functions with formulas, graphs, and mathematical tables. 10th ed. Dover; 1974.
712
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Computation of temperature distributions in transformer covers due to high crossing currents in bushing regions D. Cahue-Diaza, S. Maximova,
T
⁎,1
, R. Escarela-Perezb, J.C. Olivares-Galvanb, J. Alvarez-Ramirezc
a
Instituto Tecnologico de Morelia, PGIIE, Av. Tecnologico No. 1500, Lomas de Santiaguito, 58120 Morelia, Mich., Mexico Departamento de Energia, Universidad Autonoma Metropolitana-Azcapotzalco, 02200 Mexico City, Mexico c Departamento de Ingenieria de Procesos e Hidraulica, Universidad Autonoma Metropolitana-Iztapalapa, Mexico City 09340, Mexico b
ARTICLE INFO
ABSTRACT
Keywords: Power transformer Temperature distribution Analytical methods Heat equation Finite element method
Analytical formulae are widely applied in transformer design due to their easy implementation and low demand of computational resources. One of the great challenges in the field of transformer design is the acquisition of closed-form formulae for calculating temperature distributions in transformer tanks, taking account of different types of surrounding media. The presence of these media results in a variety of boundary conditions when solving the heat equation. In this research, a new rigorous analytical formula for calculating temperature profiles in bushing regions of transformer tank walls is obtained by solving the heat equation, with boundary conditions that properly model the heat interchange between different surrounding media. The heat source is modeled using a formula taken from a previously published investigation performed by the authors. Several study cases are performed, where analytical results are compared with finite element simulations for different boundary conditions. Analytical and numerical results show a good match, demonstrating high precision of our developed analytical method. Thus, this method provides a powerful tool for calculating temperature distributions in transformer tanks, which can be employed for optimization purposes during transformer design stages.
1. Introduction Power losses, that are originated in windings, core, tank covers, etc., result in heating of transformers. Overheating and presence of hot-spots can damage vital parts of a transformer, affecting its optimal operating regime and resulting in possible failures. Therefore, in order to improve transformer design and prevent undesirable consequences of such critical situations, it is vital to have means to compute temperature distributions and finding such hot-spots [1–5]. Several techniques for computing temperature distributions in power transformers have been developed throughout the history of transformer manufacture. The commonly used Finite Element Method (FEM) is a powerful and reliable tool for computing temperature profiles in active parts of transformers [5–9]. In particular, temperature distributions in tank walls [5,10] and oil temperature in power transformers [11–15] can be estimated by this method. In [6], FEM was employed to calculate the temperature distribution in a dry-type transformer. Joule losses in windings and eddy-current losses in magnetic cores were experimentally assessed. The results were employed to find hot-spots and air/oil flows inside the tank. The
authors emphasized the high accuracy level of the proposed experimental approach. In [7], FEM permitted the study of heat transfer and oil flow distribution, inside a 3-phase power transformer, by solving the continuity, momentum and energy balance equations in steady state. Different geometries and flow rates were analyzed in order to optimize the transformer cooling. In [8], fluid flow and convection heat transfer, in a thermo-convective circuit with several cooling channels and different inlet fluid velocities, were studied in a 3-phase power transformer immersed in mineral oil. The optimal velocity, that keeps temperature within stipulated limits, was obtained. In [9], the development of an advanced 3D finite element model, for the coupled solution of the heat transfer and fluid flow equations, that govern the transformer thermal performance, was presented. One of the main advantages of the proposed method is that there is no need to predefine convection coefficients at the interfaces between the active tank walls and circulating oil, since they are computed while the coupled heat-magnetic problem is being solved. In [12,14,15], numerical and experimental studies on a radiator-fan assembly of a power transformer were performed for different cooling configurations, to determine their cooling capacity and optimize their design.
Corresponding author. E-mail addresses: [email protected], [email protected] (S. Maximov). 1 On sabbatical leave: Universidad Autonoma Metropolitana-Azcapotzalco, 02200 Mexico City, Mexico. ⁎
https://doi.org/10.1016/j.ijepes.2019.06.019 Received 26 January 2019; Received in revised form 4 April 2019; Accepted 6 June 2019 Available online 14 June 2019 0142-0615/ © 2019 Elsevier Ltd. All rights reserved.
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Nomenclature
q q cr t kt T Ta W r r
convection coefficient global convection coefficient symmetric component of the temperature asymmetric component of the temperature Kronecker delta internal disk radius external disk radius metallic plate thickness steel electric conductivity steel relative permeability frequency
c
hc T+ T
heat flux by conduction heat flux by convection time thermal conductivity temperature ambient temperature volumetric heat source radiation emissivity Stefan-Boltzmann constant
nm
a b h
µr f
However, FEM requires considerable time-consuming computational resources, especially when dealing with 3D geometries and very pronounced skin-effects, where a high density meshing is normally required. Although numerical methods provide valuable information about thermal characteristics of transformers under fixed conditions (geometry, boundary conditions, etc.), they cannot provide information of how the temperature distribution depends explicitly on each parameter of the transformer. On the other hand, analytical formulae explicitly contain all the parameters of the transformer and surrounding media. This fact makes analytical methods powerful tools for design optimization of transformers using its parameters (geometry, materials involved in the transformer construction, etc.). Another aspect, where analytical approaches are useful, is involved with the determination of heat transfer coefficients. It is well-known that analytical or experimental determination of these coefficients is a difficult task [16,4]. Once an analytical solution to the heat equation is obtained, it can be compared with measurements. Analytical models can be fitted to experimental points by variation of heat transfer coefficients in the analytical solution, using least squares [16] or multiobjective (MO) optimization procedures [4]. Heat transfer coefficient values, thus obtained, can be consequently employed for transformer design purposes. Recently, an analytical formula to calculate temperature distribution in the bushing region of transformers has been obtained [16]. The heat source, employed to calculate the temperature profile, can be considered numerically or analytically. The proposed formula agrees well with FEM calculations. Nevertheless, it was assumed in [16] that the tank wall was surrounded by an homogeneous medium (e.g., air) with the same ambient temperature on each side and, therefore, the same heat transfer coefficient throughout the boundary between the metallic plate and the medium. Hence, the assumptions that were adopted in that work are far from real. In this paper, a new mathematical model for calculating temperature profiles in bushing regions of power transformers, that takes account of different media on each side of the tank wall, along with their heat transfer mechanisms (through boundaries), is derived. The influence of media on different sides of the wall is achieved by means of ambient temperature and heat transfer coefficients through the interface. The heat transfer coefficients depend on three mechanisms: conduction, convection and radiation. Several study cases are carried out in order to show that a good match between the new formula and FEM results is achieved for different values of the ambient temperature inside the tank. Here, the implementation of the analytical solution in the determination of the heat transfer coefficients through MO optimization procedure is discussed.
tank cover producing a nonuniform temperature distribution in the bushing region as a result of three processes: Ohm heating, thermal conduction and cooling through the surface. The tank geometry in the bushing region is presented in Fig. 1. Fig. 1B shows the cross-section of the cover of the transformer that passes along the conductor. The schematic view of this cross-section is presented in Fig. 1C. Fig. 2 shows the boundaries through which the transformer tank cover interchanges heat with different surrounding media. The upper surface 3 of the plate is in contact with air, so that the two heat transfer mechanisms through this boundary, namely convection and radiation, are to be considered. 4 is the inner surface of the transformer tank, through which the tank wall interacts with the internal medium of the transformer. The tank can be filled with air or oil, so that the thermal interaction of the tank cover with the fluid inside the transformer is given by convection heat transfer. The space inside the hole, which is a dielectric, comes into contact with the tank wall through surface 1, whereas the medium on the side of boundary 2 is air. The heat equation that governs the heat transfer phenomenon in steady state (see [16]) is given by
1 T r r r r
+
2T
z2
=
1 W r, z , kt
(1)
where the time-averaged heat source W (r , z ) can be calculated numerically [17] or analytically [18]. Due to the symmetry of the problem, cylindrical coordinates are employed. Boundary conditions on each surface can be expressed as follows:
nk ·q cr
k
= nk · q
k,
where q is the heat flux density in the tank plate, which is expressed by
2. The model A transformer bushing is an insulator that allows an electrical conductor to pass safely through the tank wall. Alternating current in the conductor produces an electromagnetic field around it, which, in turn, generates eddy currents in the tank wall. These currents heat the
Fig. 1. Bushing region of a transformer. 700
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
kt
Fourier’s law as:
kt T ,
c (T )(T
Ta) +
r
r
(T 4
Ta4 )} n ,
q cr =
(T ) +
r
T 3 mTam
r
T
T± r , z =
Ta n.
The following heat transfer coefficient is therefore obtained:
1 r r
3
(T ) +
r
T3
r
mT m a
m =0
and can be considered as a constant parameter for a given range of temperatures. Then, we find that
q cr = hc (T
Ta ) n.
kt
kt
h c1 ( T
Ta1 )
= 0, r=a
T + hc 2 (T r
Ta2 )
T + h c 3 (T z
Ta3 )
= 0, r=b
= 0, z = h /2
T (r , z ) ± T (r , 2
(r ) + T+ r
z)
1 W kt
=
T+ r
hc1 (T+
Ta1 )
kt
T+ r
+ hc 2 (T+
Ta2 )
+ h+ (T+
Ta +)
1 r r
(r ) + T r
2T
T r
hc1 T
T kt r
+ hc 2 T
kt
kt
(3)
2T + z2
kt
T kt z+
Experimental determination of the average parameter hc is discussed in Appendix B. Now, the four boundary conditions take the following form:
T kt r
(6)
.
(7)
Since the heat source W (r , z ) is an even function of z (see [18]), Eq. (1), along with boundary conditions (3)–(6), can be divided into the following two problems:
m=0
h c (T ) =
= 0, z = h /2
where T+ (r , z ) and T (r , z ) are the symmetric and antisymmetric components (with respect to the coordinate z) of the temperature distribution, which are defined in the following way:
(2)
where r is the Stephan-Boltzmann constant and n is a unitary vector normal to the surface. In [16], a linearization procedure for the convection heat flux density has been introduced. According to this method, Eq. (2) should firstly be cast in the form: 3
Ta 4 )
T (r , z ) = T+ (r , z ) + T (r , z ),
and q cr consists of two convection and radiation heat exchange terms:
q cr = {
h c 4 (T
where the existence of four different heat transfer coefficients on each surface and four different ambient temperatures near each side of the tank wall has been assumed. In order to solve Eq. (1) with the asymmetric boundary conditions (3)–(6), it is appropriate to divide the temperature field T (r , z ) into two parts as follows:
Fig. 2. Tank geometry.
q=
T z
T z
z2
+ h+T
r=a r=b
r, z , = 0, = 0,
z = h /2
=
h T |z = h /2 ,
(8)
= 0,
r=a r=b
= 0, = 0,
z = h /2
=
h (T+
Ta )z = h /2 ,
(9)
where (4)
(5)
h± =
h c3 ± h c 4 , 2
(10)
Ta ± =
hc 3 Ta3 ± hc 4 Ta 4 . h c3 ± h c 4
(11) Fig. 3. Schematic behavior of the left-hand part of (29) and solutions for qn .
701
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
where T±[N ] (r , z ) is the N th iteration for the temperature T± (r , z ) . Here, the zero iteration for the antisymmetric part of temperature distribution is taken as T [0] (r , z ) 0 . The scheme (13) and (14) can be truncated at a given iteration that provides sufficient precision of the approximate solution. The hypothesis of this work is that, with a good precision, the iteration scheme given by (13) and (14) can be considered up to the first iteration. Further, this hypothesis is verified by comparing analytical results with numerical computations.
Problems (8) and (9) are coupled through boundary conditions, so they cannot be solved separately. In the next section, the scheme of solution of these two coupled problems is developed. 3. Proposed iteration scheme Problems (8) and (9) can be decoupled if the right-hand part of the last equation of problem (8) is neglected. This term can be estimated using the mean value theorem, according to which
z
h/2, h/2 :
T r, z =
T (r , z )
T (r ,
z)
2
T = z
3.1. First iteration for the symmetric part
z.
The problem for the first iteration of the symmetric part of the temperature distribution T+ (r , z ) can be obtained from (13) by substituting N = 1 and T [0] (r , z ) 0 . The problem thus obtained displays similarity to the problem solved in [16]. However, in this case different heat transfer coefficients hc1, hc 2 and ambient temperatures Ta1 and Ta2 near the boundaries 1 and 2 are now accounted for. Following the mathematical steps given in [16], we obtain:
z=z
The following estimation can therefore be demonstrated:
h T r, z
h
T h h h= z kt
max
h /2 z h /2
max
h /2 z h /2
qz ,
(12)
where max h /2 z h /2 qz is the maximum heat flux density in the z-axis direction. Hence, this term is small if the wall thickness h is also small enough and the difference of ambient temperatures above and below the tank wall is not large. The last statement means that the heat flux density in the z-direction is small (Fourier’s law). Nevertheless, when overheating of the medium below the cover plate takes place, the right-hand part of the last equation of problem (8) becomes significant and, consequently, it cannot be neglected. In this case, an iteration scheme can be applied. In the first iteration, the righthand part of the last equation of (8) is neglected. As a result, problem (8) is solved independently of (9). In the next step, the solution of this problem is substituted into (9). In the second iteration, the solution to (9) is substituted into (8), where the right-hand part of boundary conditions is now taken into account. These consecutive steps can be represented by the following iteration scheme: 1 r r
kt kt kt
1 r r
r
[N ] T+
r
[N ] T+
2T [N ] + z2
hc1 (T+[N ]
r
1 W kt
=
= Ta + + Wn ( ) d
r
+ hc 2 (T+[N ]
z
+ h+ (T+[N ]
= 0, =
+
2T [N ]
z2
kt
T [N ] r
hc1 T [N ]
kt
T [N ] r
+ hc 2 T [N ]
T [N ] kt z
+ h+T
h T [N
1]
r , h/2 ,
(13)
= 0,
+
1 kt
b a
Gn r ,
Dn(12) I0 ( +
(11) n ) + Dn K 0 ( n (11) (22) Dn Dn Dn(12) Dn(21) (12) Dn hc 2 a2 h sin n , 2 Dn(12) Dn(21)
)
Wn ( ) d
n
2 Zn
2
(16)
Dn(11) kt
b
Dn(22) I0 (
a a1
+
(21) n ) + Dn K 0 ( n (11) (22) Dn Dn Dn(12) Dn(21) Dn(11) hc 2 a2 h sin n , 2 Dn(12) Dn(21)
Dn(11) = kt
n I1 ( n a)
Dn(12) Dn(21) Dn(22)
= kt
n K1 ( na )
= kt
n I1 ( n b)
= kt
n K1 ( nb)
)
Wn ( ) d
n
2 Zn
2
(17)
hc1 I0 ( na), + hc1 K 0 ( na),
+ hc 2 I0 ( nb), (18)
hc 2 K 0 ( nb).
Here,
= 0,
Gn (r , ) = I0 ( nr ) K 0 (
r=b
=
b
Dn(11) Dn(22)
r=a
z = h /2
nr
(15)
a
Dn(21) hc1
= 0,
[N ]
+ Bn K 0
Bn =
Ta +)
nr
,
Dn(22) hc1 a1 Dn(11) Dn(22)
z = h /2
T [N ] r
Dn(22) kt
=
r=b
[N ] T+
An I0
An
= 0,
Ta2 )
nz
where
r, z ,
Ta1 )
cos n=0
r=a
[N ] T+
r
+
T+[1] r , z
h
(T+[N ] (r ,
h/2)
Ta ),
n
) (
r ) + I0 (
n
) K 0 ( nr )[1
(
r )] (19)
(14)
is the Green’s function, where I0 ( nr ) and K 0 ( nr ) are modified Bessel
Fig. 4. Mesh used for FEM analysis. 702
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Finally, {
Table 1 Parameters for simulations. Parameters
Case I
Case II
Case III
µr kt f h c1 h c2 h c3 h c4 Ta1 Ta2 Ta3 Ta4
500 16 W/mK 50 s−1 25 W/m2K 10 W/W/m2K 10 W/W/m2K 35 W/W/m2K 80 °C 18 °C 18 °C 80 °C
500 16 W/mK 50 s−1 25 W/W/m2K 10 W/W/m2K 10 W/W/m2K 35 W/W/m2K 100 °C 18 °C 18 °C 100 °C
500 52 W/mK 50 s−1 25 W/W/m2K 10 W/W/m2K 10 W/W/m2K 35 W/W/m2K 100 °C 18 °C 18 °C 100 °C
n
{0, }} are solutions to the following equation:
n
h h = +. 2 kt
tan
(21)
Eq. (21) can be solved either numerically or analytically (see [16]). Here, the asymptotic analytical solution to Eq. (21), obtained in [16], is employed. 3.2. First iteration for the antisymmetric component Solution (15) should be substituted into the right-hand side of the last equation of (9) with N = 1. The solution to (9) can be sought using separation of variables [16] as follows: (22)
T [1] (r , z ) = R (r ) Z (z ). As a result, the following equations can be obtained:
R +
1 R + q2R = 0, r
(23)
Z
q2 Z = 0
(24)
with boundary conditions
[kt R hc1 R]r = a = 0, [kt R + hc 2 R]r = b = 0.
(25)
General solutions to Eqs. (23) and (24) can be obtained in the form:
Wn (r ) =
2 h + sin( nh)/
h /2 n
h /2
nz
dz.
Z = sinh(qz ),
(27)
(kt qJ1 (qa) + hc1 J0 (qa)) C1 + (kt qY1 (qa) + hc1 Y0 (qa)) C2 = 0,
r ) is Heaviside unit-step function.
W r , z cos
(26)
where J0 (qr ) and Y0 (qr ) are Bessel functions of zero order and q is an unknown constant. Substitution of the general solution (26) into the boundary conditions (25) yields the following linear system of algebraic equations:
Fig. 5. 2D temperature distribution. FEM-simulations. Case I.
functions of zero order, while (
R = C1 J0 (qr ) + C2 Y0 (qr ),
(kt qJ1 (qb)
hc 2 J0 (qb)) C1 + (kt qY1 (qb)
hc 2 Y0 (qb)) C2 = 0.
(28)
The existence of a nontrivial solution to system (28) is provided if the determinant of the coefficient matrix of system (28) is equal to zero:
(20)
Fig. 6. Temperature profiles on the upper surface: theoretical and FEM-calculated curves. Case I. 703
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Fig. 7. Temperature profiles on the lower surface: theoretical and FEM-calculated curves. Case I. n (r )
Det = (kt qJ1 (qa) + hc1 J0 (qa))(kt qY1 (qb) (kt qJ1 (qb)
hc 2 Y0 (qb))
hc 2 J0 (qb))(kt qY1 (qa) + hc1 Y0 (qa)) = 0,
= (kt qn Y1 (qn a) + hc1 Y0 (qn a)) J0 (qn r ) (29)
} are orthogonal, as shown in Appendix A. The functions { n (r ) n The norm of these functions is also calculated there. Therefore, instead of this system, the following set of orthogonal and unitary functions can be considered:
which is the base to find the values of constant q. Eq. (29) has an in} , which can be numerically obtained. finite number of roots {qn n The solutions for q are shown in the Fig. 3. After substituting qn into (28) and taking account of (29), the following general solution to system (28) is obtained:
C1, n = Cn (kt qn Y1 (qn a) + hc1 Y0 (qn a)) C2, n = Cn (kt qn J1 (qn a) + hc1 J0 (qn a)).
n (r )
n (r )
=
n (r )
.
(32)
n
As a result, the general solution for the asymmetric part of the temperature distribution takes the form:
(30)
Thus, the solutions to Eqs. (23) and (24) take the form:
Rn (r ) = Cn
(kt qn J1 (qn a) + hc1 J0 (qn a)) Y0 (31)
(qn r ).
T [1] r , z =
and Zn (z ) = sinh(qn z ),
Cn
n (r )sinh
qn z .
n=1
(33)
Substitution of this solution into the last boundary condition of problem
where
Fig. 8. Temperature profiles on the upper and lower surfaces. Analytical solution. Case I. 704
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
T [1] r , z h sinh(qn z )
= T+[1] r , z + n= 1
×
(9) yields:
Cn
n (r )(qn kt cosh(qn h /2)
=
h (T+[1] (r , h/2)
+ h+sinh(qn h/2))
h
b a
(Ta
T+[1] ( , h /2))
n(
) d
qn kt cosh(qn h/2) + h+sinh(qn h /2)
.
T+[1] ( , h/2))
n(
) d .
(36)
Three different study cases are considered in this section to illustrate the accuracy of our analytical approach. Ansys Fluent software has been employed to carry out FEM simulations considering an axisymmetric geometry (see Fig. 2). A cylindrical metallic plate with internal and external radii a = 0.0125 m and b = 0.265 m, respectively, and a thickness h = 0.01 m, has been considered. These geometry parameters were used for all the simulations. The cross-section of the tank wall was covered using a mesh with a refinement near the plate surface, as shown in Fig. 4, so as to accurately compute eddy current losses (the power density of the eddy current losses is higher near the tank wall surface [18,19]). The electromagnetic model reported in [18] was employed for calculating eddycurrent losses while our new analytical method is used to compute the temperature distribution. Since the model [18] faithfully reproduces the losses, there was no need to consider different currents in the
Using the orthogonality property of the system of functions { n (r ) n } , the following expression for the unknown constants Cn can be obtained from Eq. (34):
Cn =
(Ta
4. Study cases (34)
Ta ).
a
Eq. (36), together with (15), represents a new model for calculating temperature distributions in bushing areas of transformer tanks. The implementation of this model requires: the heat source function W (r , z ) (calculated numerically or analytically), the geometry parameters, properties of the tank wall material and surrounding fluid and the heat transfer coefficients on the tank wall surface. The algorithm for calculating the temperature distribution consists of the following steps: (i) calculation of average ambient temperatures (11) and heat transfer coefficient (10) on surfaces 3 and 4 ; (ii) estimation of the parameters n by solving Eq. (21) numerically or analytically (see [16]); (iii) computation of Fourier coefficients (20) of the heat source; (iv) calculation of the parameters (18) and constants (16) and (17); and (v) finally, implementation of formula (36). Finally, the appropriate number of terms that should be taken from series (36), to get a good precision, is estimated in Appendix C.
Fig. 9. 2D temperature distribution. FEM-simulations. Case II.
n=1
n (r )
b
qn kt cosh(qn h /2) + h+sinh(qn h/2)
(35)
Substitution of (35) into (33) and summing the symmetric (15) and antisymmetric (33) components of the solution, the following temperature distribution in the first order iteration is obtained:
Fig. 10. Temperature profiles on the upper surface: theoretical and FEM-calculated curves. Case II. 705
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Fig. 11. Temperature profiles on the lower surface: theoretical and FEM-calculated curves. Case II.
conductor to represent different levels of power loss. Therefore, all the simulations were performed using the same value I = 1000 A of the current in the conductor, i.e., the same volumetric heat source was considered in each case. So, the simulations of our new work were focused on modeling different boundary conditions round the tank wall, namely, different ambient temperatures. The operating conditions and parameter values of the system are given in Table 1. The loss density employed for the analytical computations of the temperature distributions was calculated using the electromagnetic model published by the authors in [18]. Numerical computations were performed using FEM both for calculating power losses and for computing the temperature distribution. The results of the simulations are presented in Figs. 5–16. Case I corresponds to the temperature of the medium below the transformer cover Ta4 = 80 °C and thermal conductivity of the tank wall kt = 16 W/mK . Fig. 5 shows the axisymmetric temperature distribution obtained on the upper surface of the cover, using FEM-simulations. Figs. 6 and 7 compare the results of FEM-simulations with our
analytical results, obtained with formula (36), on the upper and lower surfaces of the tank wall, respectively. It can be observed in Figs. 6 and 7 that both numerical and analytical results closely match. In Fig. 8, the difference between analytically calculated temperatures on the upper and lower surfaces, due to difference between ambient temperatures above and below the tank wall, can be observed. In the second case, the temperatures below the covering plate and in the hole were set to Ta4 = 100 °C and Ta1 = 100 °C , in order to model oil overheating in the transformer and providing a higher temperature gradient in the tank wall thickness. The numerical and analytical results are presented in Figs. 9–12. It can be observed in these figures that the increment of 20 °C in the temperature near the lower surface leads to increased temperature values at each point of the tank wall. However, it did not noticeably affect the temperature gradient in the tank wall thickness. The third case (Figs. 13–16) showed that the increment of the thermal conductivity of the tank wall from 16 W/mK to 52 W/mK considerably decreases the temperature gradient in both the axial and
Fig. 12. Temperature profiles on the upper and lower surfaces. Analytical solution. Case II. 706
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
exchange through surface 2 is not an important factor, since the area of this surface is small with respect to 3 and 4 (the plate thickness h is small with respect to the tank dimensions). Moreover, the heat-exchange velocity near 2 is much smaller than in the bushing region, since the temperature near the outer border is smaller. Therefore, the surface 2 geometry can be disregarded. However, the total area of boundaries 3 and 4 is an important factor in the heat transfer problem. In this section, it is shown that our analytical development can be generalized for actual geometries such as that presented in Fig. 1, where the bushing is displaced from the center of the plate. It is assumed that the radius of the area, where the main part of losses is concentrated (e.g., 90% of total losses), is much smaller than the outer dimensions of the cover. To do so, we consider formulae (15)–(19). In virtue of the asymptotic behavior of the modified Bessel functions (Kn ( nr ) rapidly decreases whereas In ( nr ) increases as the distance r grows [20]), it follows from (18) that 1/ Dn(21) and Dn(22) are small parameters. Formally, 0 . Application of this fact to Eqs. it means that Dn(21) and Dn(22) (16) and (17) shows that the parameter An is small enough to be neglected in Eq. (15) and the parameter Bn can be simplified as follows: Fig. 13. 2D temperature distribution. FEM-simulations. Case III.
Bn
radial directions. The three analyzed cases demonstrate that our analytical formula (36) is reliable for calculating temperature distributions and evaluating hot spots in bushing regions of transformer tank walls. The selected iteration scheme (13) and (14), truncated at the first order iteration resulted in a good precision when computing temperature distribution even in the case of an extremely high temperature below the tank wall.
1 Dn(11) kt Dn(12)
b a
K0
Wn ( ) d +
n
n
2 Zn
2
hc1 a1 sin Dn(12)
nh . 2
(37) As a result, solution (15) takes the form:
T+[1] r , z = Ta + +
cos
nz
{
n=0
5. Generalization for nonaxisymmetric geometries
+
Results (15) and (36) were obtained for a geometry with axial symmetry. However, bushings are mounted at points displaced from the center of the tank cover in many cases. Therefore, applicability of the obtained results for real cases, such as that shown in Fig. 1A, should be verified. Otherwise, solution (36) should be generalized for other geometries, the path that is undertaken next. It is well-known that the heat source is mainly concentrated in an area near the conductor, so that the main contribution to the tank wall heating is localized in the bushing region. It is also clear that the heat
n
2 Zn
2
h c1 a 1 sin Dn(12)
1 Dn(11) kt Dn(12) nh
2
K0
b a
K0
nr
Wn ( ) d
n
+
1 kt
b a
Gn r ,
Wn ( ) d }. (38)
The analysis of this formula shows that: (a) formula (38) does not depend on the ambient temperature T2 , on the right side of the covering plate; (b) the dependence of (38) on the outer shape of the transformer cover is only given by the upper limit of the integrals in the right-hand side of (38). Hence, the dependence of (38) on the covering plate shape is very weak. It can be proved by taking a derivative of (38) with
Fig. 14. Temperature profiles on the upper surface: theoretical and FEM-calculated curves. Case III. 707
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
Fig. 15. Temperature profiles on the lower surface: theoretical and FEM-calculated curves. Case III.
respect to the outer radius b that
T+[1] = b
cos n=0
nz
{
1 Dn(11) K0 kt Dn(12)
where Eq. (40) represents a convex curve. Thus, solution (41) can be generalized as follows: nr
K0
nb
+
1 Gn r , b } Wn (b) b, kt
T+[1] r , z
(39) which is negligible, since the heat source W (r , z ) (and consequently the Fourier coefficients Wn (r ) ) decreases rapidly by increasing the distance r. Moreover, the values of the Bessel function K 0 ( nb) on the cover border b is also negligible. This and the fact that Eq. (38) does not depend on the temperature near the boundary 2 and the convection coefficient hc2 , proves that solution (38) depends weakly on the boundary 2 shape. In turn, this allows the generalization of solution (38) for the case of a transformer cover of an arbitrary geometry. Let us consider a transformer cover with a shape given by the following equation (see Fig. 17):
r = b ( ),
[0, 2 ],
= Ta + +
cos
nz
{
n=0
+ n b( ) a
2 Zn
2
h c1 a 1 sin Dn(12)
Gn r ,
Dn(11) 1 Dn(12) 2 kt nh 2
K0
2 0
nr
Wn ( ) d },
d
+
b( ) a
1 2 kt
K0 2 0
n
Wn ( ) d
d
(41)
where the integrals in the right-hand side are integrated over the boundaries 3 and 4 . In the case of b = const , Eq. (41) is reduced to (38).
(40)
Fig. 16. Temperature profiles on the upper and lower surfaces. Analytical solution. Case III. 708
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
result, we obtain:
T [1] r , , z =
T [N ] 1 1 r + 2 r r r r
2T [1] 2
+
2T [N ]
z2
T [N ] 1 r + r r r
2T [N ]
z2
= 0,
which proves that generalization (42) satisfies the heat equation even for the nonaxisymmetric case. In conclusion, it should be noticed that b = const means that Eq. (42) takes the form (36). Thus, a good precision solution, generalized for the case of an arbitrary geometry, takes the form:
T [1] (r , , z ) = T+[1] (r , , z ) + T [1] (r , , z ), where T+[1] (r , , z ) and T [1] (r , , z ) are represented by (41) and (42), respectively. 6. Conclusion Fig. 17. Asymmetric geometry of the transformer cover.
In this paper, a new analytical formula for calculating temperature distributions in bushing regions of transformer tanks has been obtained. The new formula takes account of the fact that the tank wall could be surrounded by different media with different thermal properties such as the heat transfer coefficient, ambient temperature and heat radiation emissivity. In particular, the plate can be in contact with cooling fluids on the inner side (e.g., oil) and air on the outer side. Unlike previous works, all these possibilities have been taken into account in our proposed formula, through lumped parameters of heat transfer coefficients and different ambient temperatures. As a result, the new analytical approach is capable to reproduce faithfully temperature gradients in the tank wall thickness, produced by the difference of ambient temperatures above and below the tank wall, which was impossible with previous analytical models (see, e.g., [16,4]). The analytical results have been compared with FEM simulations, which have also been performed considering heat source calculations based on electromagnetic FEM calculations. The comparison shows that our new formula for temperature distributions leads to a good match with numerical results. In turn, this demonstrates that the proposed analytical method is reliable for predicting temperature distributions in transformer tanks and estimating hot spots. The analytical formula provides high precision results with no need of expensive computational resources. Unlike numerical methods, the analytical formula explicitly depends on all the parameters of the system such as geometry and material parameters involved in the transformer construction. This makes our analytical method preferable for transformer optimal design with respect to its parameters and for determination of thermal parameters through MO optimization procedures.
A similar arbitrary geometry generalization for the antisymmetric part of (36) can be realized as follows:
T [1] r , , z
n=1
h sinh(qn ( ) z ) × qn ( ) kt cosh(qn ( ) h/2) + h+sinh(qn ( ) h/2)
1 2
0
=
2
d
b( ) a
(Ta
T+[1] ( , h/2))
,
n
n
r,
d .
(42)
Here, the dependence of the parameters qn on the angle can be obtained from (29), where parameter b is a function of (see (40)). Finally, functions n (r , ) , as well as Eq. (42), depend on the angle through parameters qn and b. Let us show that (42) satisfies the heat equation with good precision. To do so, the derivative of (42) with respect to should be estimated:
T [1]
T [1] qn T [1] b + = qn b
=
T [1] dqn T [1] n db + . qn db d n db
(43)
In order to estimate the derivative of dqn / db , a derivative of Eq. (29) with respect to b should be taken. As a result, after some bulky, but easy calculations, the following estimation can be obtained:
dqn db
qn
=
+
b
1 , b2
(44)
where is the Landau big O notation. In other words, parameters qn depend weakly on the cover shape. In a similar way, the estimation of the derivative of n (r ) with respect to b shows that it is negligible. Summarizing, we are led to the conclusion that derivative (43) is small
(1/b 2)
enough. Similarly, it can be shown that
2T [1] 2
Acknowledgements The authors are grateful for financial support provided by the following Mexican Institutions: CONACYT, SNI and PRODEP.
is also negligible. As a
Appendix A
} can be proved in a standard way. Let us consider two functions of this set: The orthogonality of the set of functions { n (r ) n These functions satisfy Eq. (23) for different values of q and boundary conditions (25): m
+
1 r
m
n
+
1 r
n
where qm
+ qm2
+ qn2
m
n
= 0,
and
n (r ) .
(45)
= 0,
qn . Multiplication of Eq. (45) by r
m (r )
(46) n (r )
and Eq. (46) by r
m (r )
709
and subtraction of the second equation from the first one, leads, after
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
simplifying and integrating with respect to r over the segment [a , b], to the following equation:
r(
n
m
b n )|a
m
+ (qm2
b
qn2 )
a
m
n rdr .
(47)
Substitution of boundary conditions (25) into (47) leads to the following orthogonality condition b
n rdr
m
a
=0
n. for m In order to find the norm of functions (q, r ) = (kt qn Y1 (qn a) + hc1 Y0 (qn a)) J0 (qr ) instead of m , where q can be obtained: n
b
2 n rdr
a
=
r(
lim
q
n q2
qn
qn2
let us substitute into Eq. (47) the function (48)
(kt qn J1 (qn a) + hc1 J0 (qn a)) Y0 (qr )
qn . It is evident that
2
=
n (r ) ,
b n )|a
=
(q, r ) satisfies Eq. (23) and
1 lim r 2qn q qn
(qn , r ) =
n (r ) .
Then, the following expression for the squared norm
b
q
n
q
,
n
(49)
a
where L’Hôpital’s rule has been used to calculate the limit. From (48), it can be directly verified that
q
=
qr
and
=
q
r q
.
(50)
Substitution of the properties (50) into (49) yields: 2
=
n
r2 (q 2 2qn2 n
2 n
+
2 b n )|a .
This result, with boundary conditions (25) taken into account, finally takes the form: 2
=
n
h2 1 b2 qn2 + c22 2 2qn kt
2 n (b )
a2 qn2 +
hc21 kt2
2 n (a )
.
(51)
Appendix B The difficulty of the exact theoretical determination of thermal parameters such as the heat transfer coefficient hc , radiation emissivity, etc., is well known (see [4,5,16]). Moreover, convection coefficients are, in general, temperature-dependent functions. Also, a direct determination by experiment is often complicated [21]. However, temperature assessment of transformers is sensitive to the precision of these parameters. Since our new formula (36) showed a good match with FE results, it can be effectively used for calculating thermal parameters. Recently, multi-objective optimization was efficiently employed for thermal parameter determination [4] through comparison of FEM simulations with measurements. The precision of our new formula permits its implementation within MO optimization, where thermal parameters hc1, hc 2, hc3 and hc4 are employed as independent variables. Let us denote Tk as the temperature values measured at different points on the surface of the metallic disc, with coordinates rk , where 4 is the following vector: k = 1, 2, …, N , by N sensors, whereas analytical results at the same points, obtained from (36), are: T (rk, H) , where H
H = (hc1, hc 2 , hc3, hc 4 ). Thus, the first objective function to be minimized is established as the N
f1 (H) =
k=1
2 -space
error:
Tk ) 2 .
(T (rk, H)
(52)
The minimization of only one parameter (52) with respect to the vector H (the so-called single-objective optimization) is the widely-known leastsquare method. Another objective function that was proposed in [4] is the sensitive function that is determined as follows:
f2 (H) =
max1
i 4 f1i
f1
f1
,
(53)
where f1i is the first objective function (52) corresponding to ± 4% relative perturbations of the parameters hci , where i = 1, 4 . Thus, the 2D vector f = (f1 , f2 ) should be optimized with respect to the 4D vector H . The advantage of the use of the analytical solution instead of FEM simulations in determining thermal parameters is that the analytical calculations are precise and do not require high-consuming computational resources. Appendix C The right-hand side of (36) is a generalized Fourier series with orthogonal and unitary system of functions { Parseval’s identity, namely, by the convergence of the series
n (r )} .
Its convergence is provided by
2
cn
<
,
(54)
n= 1
where 710
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al.
cn =
h sinh(qn z ) qn kt cosh(qn h/2) + h+sinh(qn h/2)
(55)
are coefficients of the series in (36). The number of terms that should be taken into account in (36) depends on the required precision, and it can be estimated based on series (54). Firstly, the minimum number of terms N is so that qN is big enough. This fact allows the following estimation for coefficients (55):
h sinh(qn h /2)
cn
qn kt cosh(qn h/2) + h+sinh(qn h/2)
=
h qn kt coth(qn h/2) + h+
h . qn kt
(56)
Therefore, series (54) is majorized by another one as follows: 2
cn n= 1
n= 1
h2 qn2 kt2
(57)
Secondly, parameters qn can be asymptotically estimated from (29) for the high values of n. Substitution of the asymptotic formulae for Bessel functions [20] into Eq. (29), after neglecting small terms, yields:
sin(qn (b
a))
0.
The solutions to this equation are:
n
qn
b
a
,
where n 2
cn n= 1
(58) . Thus, the convergence of series (57) can be estimated in the following way:
h2 kt2
n=1
1 qn2
h 2 (b kt2
a) 2 2
n=1
1 , n2
where
n= 1
2 1 = . 2 n 6
(59)
Let us truncate series (59) at the
rN = n=N +1
1 = n2
N th
term. Then, the residual term rN of the truncated series is:
N+1 ,
(60)
where
N+1 =
d ln (x ) dx
x=N+1
is the digamma function [22]. Therefore, in order to provide a good precision, we can make the residual term do not exceed a small part of the 1. Therefore the minimum number N can be found from the following complete series value. Namely, let be a small parameter such that 0 < equation:
N+1 For example, if N 11.
2
6
.
= 0.01 (i.e., the relative error does not exceed 1% ), the approximate value of N is 60. In the case of
= 0.05 (5% of the relative error)
[8] Chereches NC, Chereches M, Miron L, Hudisteanu S. Numerical study of cooling solutions inside a power transformer. Energy Procedia 2017;112:314–21. [9] Tsili M, Amoiralis EI, Kladas AG, Souflaris AT. 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. [10] Sitar R, Šulc I, Janić Ž. Prediction of local temperature rise in power transformer tank by fem. Procedia Eng 2017;202:231–9. [11] Najar S, Tissier JF, Cauet S, Etien E. Improving thermal model for oil temperature estimation in power distribution transformers. Appl Therm Eng 2017;119:73–8. [12] Anishek S, Sony R, Kumar JJ, Kamath PM. Performance analysis and optimisation of an oil natural air natural power transformer radiator. Procedia Technol 2016;24:428–35. [13] Alegi GL, Black WZ. Real-time thermal model for an oil-immersed, forced-air cooled transformer. IEEE Trans Power Deliv 1990;5(2):991–9. [14] Paramane SB, Van Der Veken W, Sharma A. A coupled internal-external flow and conjugate heat transfer simulations and experiments on radiators of a transformer. Appl Therm Eng 2016;103:961–70. [15] Garelli L, Ríos Rodriguez G, Storti M, Granata D, Amadei M, Rossetti M. Reduced model for the thermo-fluid dynamic analysis of a power transformer radiator working in ONAF mode. Appl Therm Eng 2017;124:855–64. [16] Maximov S, Escarela-Perez R, Olivares-Galvan JC, Guzman J, Campero-Littlewood
References [1] Allahbakhshi M, Akbari M. Heat analysis of the power transformer bushings using the finite element method. Appl Therm Eng 2016;100:714–20. [2] Cinar MA, Alboyaci B, Sengul M. Comparison of power loss and magnetic flux distribution in octagonal wound transformer core configurations. J Electr Eng Technol 2014:1290–5. [3] Díaz-Chacón JM, Hernandez C, Arjona MA. Finite element and neural network approach for positioning a magnetic shunt on the tank wall of a transformer. IET Electr Power Appl 2016;10:827–33. [4] Penabad Durán P, Di Barba P, Lopez-Fernandez X, Turowski J. Electromagnetic and thermal parameter identification method for best prediction of temperature distribution on transformer tank covers. Int J Comput Math Electr Electron Eng 2015;34:485–95. [5] Penabad-Duran P, Lopez-Fernandez XM, Turowski J. 3D non-linear magnetothermal behavior on transformer covers. Electr Power Syst Res 2015;121:333–40. [6] Ortiz C, Skorek AW, Lavoie M, Bénard P. Parallel CFD analysis of conjugate heat transfer in a dry-type transformer. IEEE Trans Ind Appl 2009;45:1530–4. [7] Wakil NE, Chereches NC, Padet J. Numerical study of heat transfer and fluid flow in a power transformer. Int J Therm Sci 2006;45(6):615–26.
711
Electrical Power and Energy Systems 113 (2019) 699–712
D. Cahue-Diaz, et al. E. New analytical formulas for temperature assessment on transformer tanks. IEEE Trans Power Deliv 2016;31(3):1122–31. [17] Olivares-Galvan JC, Perez RE, Kulkarni SV, Leon FD, Venegas-Vega MA. 2d finiteelement determination of tank wall losses in pad-mounted transformers. Electr Power Syst Res 2004;71:179–85. [18] Maximov S, Olivares-Galvan JC, Escarela-Perez R, Magdaleno-Adame S, CamperoLittlewood E. New analytical formulae for electromagnetic field and eddy current losses in bushing regions of transformers. IEEE Trans Magn 2015;51. 6 300 710–1–10.
[19] Maximov S, Olivares-Galvan JC, Magdaleno-Adame S, Escarela-Perez R, CamperoLittlewood E. Calculation of nonlinear electromagnetic fields in the steel wall vicinity of transformer bushings. IEEE Trans Magn 2015;51. 6 300 811–1–11. [20] Olver FW. Asymptotics and special functions. 3rd ed. Academic Press; 1974. [21] Kitak P, Glotic A, Ticar I. Real-time thermal model for an oil-immersed, forced-air cooled transformer. IEEE Trans Magn 2015;50(2):933–6. [22] Abramowitz M, Stegun IA. Handbook of mathematical functions with formulas, graphs, and mathematical tables. 10th ed. Dover; 1974.
712