Energy Conversion and Management 52 (2011) 2565–2574
Contents lists available at ScienceDirect
Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman
Analytical calculation of detailed model parameters of cast resin dry-type transformers M. Eslamian, B. Vahidi ⇑, S.H. Hosseinian Department of Electrical Engineering, Amirkabir University of Technology, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 12 September 2009 Received in revised form 31 July 2010 Accepted 15 January 2011 Available online 24 March 2011 Keywords: Cast resin dry-type transformer Detailed model Parameter calculation FRA FEM Impulse test
a b s t r a c t Non-flammable characteristic of cast resin dry-type transformers make them suitable for different kind of usages. This paper presents an analytical method of how to obtain parameters of detailed model of these transformers. The calculated parameters are compared and verified with the corresponding FEM results and if it was necessary, correction factors are introduced for modification of the analytical solutions. Transient voltages under full and chopped test impulses are calculated using the obtained detailed model. In order to validate the model, a setup was constructed for testing on high-voltage winding of cast resin dry-type transformer. The simulation results were compared with the experimental data measured from FRA and impulse tests. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Transformers in operation are subject to various kinds of overvoltages caused by lightning strikes, switching operations or system disturbances. For designing the insulation of a transformer suitable for all kinds of overvoltages, the voltage stresses within the windings need to be determined. For this purpose, voltage distributions within the transformer windings for the specific test voltages are calculated. Generally, under normal conditions in a steady state, the voltage is distributed linearly inside the winding. For impulse voltage, the voltage distribution in the winding is not linear and for the calculation of impulse voltage distribution in the windings, they are required to be simulated in terms of an equivalent circuit consisting of lumped R, L, C elements. In this paper, the transient model of cast resin dry-type transformers is presented. Cast resin dry-type transformers are the most suitable transformers for distribution of electricity in high degree of safety. Dry-type transformers compared with oil-immersed are lighter and non-flammable. They also do not have contaminating substances such as oil. Non-flammable characteristic of these transformers make them suitable for residential and hospital usages. Fig. 1 shows the structure of one phase of a cast resin dry-type transformer consisting of low and high-voltage windings. Low⇑ Corresponding author. Tel.: +98 21 64543330; fax: +98 21 66406469. E-mail addresses:
[email protected] (M. Eslamian),
[email protected] (B. Vahidi),
[email protected] (S.H. Hosseinian). 0196-8904/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enconman.2011.01.011
voltage winding of this transformer is made of a full-height aluminum sheet wounded simultaneously with insulation layer. Highvoltage winding is constructed from several series disks that every disk is made of aluminum foils interleaved with insulating layers. After the high-voltage winding is wounded, it is placed in a mold and cast in a resin under vacuum pressure. Lower sound levels are realized as the winding is encased in solid insulation. Filling the winding with resin under vacuum pressure eliminates voids that can cause corona. With a solid insulation system, the winding has superior mechanical and short-circuit strength and is impervious to moisture and contaminants. Transformer transient modeling has been a subject of investigation and research for a century. In the 1954, Lewis [1] proposed that the transient behavior of a transformer winding can be studied with an equivalent ladder-type network composed of a finite number of uniform sections. Lewis’s model was applicable only to a uniform winding. Furthermore, the representation of inductive coupling effect was included by modifying the self-inductance value. Subsequently with the advent of digital computers, McWhirter et al. [2], studied the same problem based on an equivalent circuit approach. Their model still suffers the restrictions arising from the size and symmetry of the equivalent circuit used to represent the winding. Dent et al. [3] used an equivalent circuit of the same general form as the proposed by Lewis, but with certain differences. Dent’s model can represent a non-uniform winding and the effect of inductive coupling between sections is taken into account. After Dent’s paper was published, most of the researchers in this area, concentrated their attention on the calculation of parameters of
2566
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
Core leg
HV disks
air duct between HV and LV windings
casted resin LV layers Fig. 1. Schematic view of cast resin transformer.
the equivalent circuits. Okuyama [4] calculated the self and mutual inductances of transformer winding through introduction of some correction factors obtained from experiment. Stein [5] and Kawaguchi [6] proposed a method to calculate equivalent series capacitance by computing the electrostatic energy stored in the coils. An equivalent network for a multi-winding transformer has been reported in [7] in which the conventional ladder network used for a single winding consisting of lumped elements is extended for multiple windings. De Leon and Semlyen [8] used the image method to calculate turn to turn leakage inductance and the charge simulation method to calculate the capacitance between turns and from turn to ground. They also proposed a detailed model of losses for accurate calculations [9]. Such a model, although very accurate, may be prohibitive from the point of time and memory of computers. Until the introduction of the computer, there was a lack a practical analytical formulas to compute the mutual and selfinductances of coils with an iron core. Rabins [10] developed an expression to calculate mutual and self-inductances for a coil on an iron core based on the assumption of a round core leg and infinite core yokes, both of infinite permeability. Fergestad and Heriksen [11] improved Rabins’s inductance model in 1974 by assuming an infinite permeable core except for the core leg. They also introduced a method of calculating voltage oscillations in multi-winding transformers during impulse test [12]. White [13] derived an expression to calculate the mutual and self-inductances in the presence of an iron core with finite permeability under the assumption of an infinitely long iron core. Wilcox et al. [14,15] derived a frequency-dependent impedance formula to calculate self and mutual impedances of coil sections on the assumption of a solid homogeneous core. Application of this formula accounts for the effects of transient flux penetration into the transformer core, including the effects of frequency-dependent inductances and transient eddy currents in the core laminations. Nowadays, numerical methods especially finite element method (FEM) has become more popular between engineers. The most important advantage of FEM is that any complicated geometry is solvable by this method because the formulation of FEM is independent of the geometry. Until now many works about important phenomena of transformers such as inductance calculation and electric field analysis have been done successfully by this method [16,17]. 3D FEM computation of the high frequency parameters using a 3D complex permeability is presented in [18,19]. Frequency-dependent losses due to eddy currents are represented in [20] by means of a frequency-dependent complex permeability for iron core. Although FEM is an accurate method but it does not manifest basic concepts of the problem. Instead, using analytical methods makes it possible to directly involve with solving of equations which will provide a deep perception of the problem. On the other hand the accuracy of the FEM solution is usually a function of the mesh resolution while it requires substantial
amounts of computer and user time. So analytical methods are still used and privileged by many designers especially where a quick theoretical solution of the problem is needed. However it will be so beneficial if the analytical results could be verified by the FEM and then modified if it is required. This is the main procedure that is used in this paper for calculation of detailed model parameters of cast resin dry-type transformers. Up to now, all works done on transformer transients are for oil filled transformers and this was not done for dry-type transformers yet. As it is clear from its name, dry-type transformer do not have any oil and instead a solid dielectric material such as resin is used in its structure. The shape of conductor, types of dielectric materials and the structure of windings of dry-type transformers are different from oil filled transformers. Also the shape of electrical field is rather different. On the other hand, the surrounding air as an insulator and cooler is running in the space between windings and tips of bars. Because of large difference of air dielectric coefficient with other solid insulators of dry-type transformers (about four times), the large amount of voltage in this kind of transformers is dropped on air spaces and with respect to low resistance of air in compare with other solid materials, it’s important to have information about the size and the way of electrical field distribution in transients to make sure about tolerance of air spaces. Therefore, study and research on dry-type transformers seems to be important and useful. In this paper we tried to study and simulate the transient behavior of dry-type transformers as a new generation of distribution transformers by using analytical formulas. Also all calculated parameters are compared with the FEM results and if it was necessary, correction factors are introduced for modification of the analytical solutions. 2. The equivalent circuit of dry-type transformer Equivalent circuit used for modeling of dry-type transformers is shown in Fig. 2. In this case, each section of LV winding which is placed between two air canals and every disk of HV winding is considered as a branch in the equivalent circuit. Every branch is composed of parallel connection of self-inductance of Li and series capacitance of Ki. There is a mutual inductance between every pair of branches (Lij, Lik, . . .). Regarding that the last layer of LV is creating an equipotential surface in front of other layers of LV, only the capacitances of the last layer of LV to HV disks are considered which are shown with Cp. The Ce is the capacitance of the first layer of LV to the core which is grounded. Ri is ohmic losses of each branch and Gi, Gp and Ge are conductance which is representing the dielectric losses. HV1
Gi Ge
Lij Ui-1 Ri
Ce LV1
Ki
L i Ui LV2
Gp
Cp HV2 Fig. 2. Equivalent lumped parameter circuit.
2567
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
3. Inductance calculation Traditionally, the magnetic flux is composed of two components of leakage and common fluxes. The leakage flux is not very dependent on frequency however the common flux which is passing through the core is highly frequency-dependent and in very high frequencies, it is completely displaced from the middle of the core and only flows on the surface of the core. Since common flux contributes to the same voltage induced in all turns and because of the reason that the most of the magnetic energy of system results from leakage flux, it is therefore the leakage inductance that is the determining factor for internal transients in transformer windings and the effect of iron core could be neglected in transformer high frequency modeling. It is especially true when the LV winding is short-circuited and end of HV winding is grounded as in our study. However, the existence of core leg change the shape of leakage flux and as a hollow, thin metallic cylinder, it will create a flux-normal boundary condition for leakage flux lines. 3.1. Calculation of inductances considering core leg effect White derived an expression for the mutual and self-inductance between sections of a transformer winding by solving a twodimensional problem in cylindrical coordinates for the magnetic vector potential assuming a non-conductive and infinitely long open core [13]. He assumed that the leakage inductance of an open-core configuration is the same as the closed core. Starting from the definition of the magnetic vector potential ~ B ¼ r~ A ~ ~ ~ and r A ¼ 0 and using Ampere’s law J ¼ r H, White solved the following equation:
Aðr; zÞ ¼ l0~ Jðr; zÞ r 2~
ð1Þ
The solution broke into two parts: the air-core solution and the change in the solution due to the insertion of the iron core, as shown in the following equations [13]:
( ~ Aðr; zÞ ¼
~ A0 ðr; zÞ þ ~ A1 ðr; zÞ; 0 6 r < Rc ~ A2 ðr; zÞ; r > Rc A0 ðr; zÞ þ ~
ð2Þ
~ A1 ðr; zÞ and A0 ðr; zÞ is the solution when the core is present and ~ ~ A2 ðr; zÞ are the solutions when the iron core is added. Applying Fourier series to Eq. (1), the solution for ~ A0 ðr; zÞ was found first and then ~ A2 ðr; zÞ. Knowing the magnetic vector potential allows the flux linking a filamentary turn at (r, z) to be determined by recalling H /ðr; zÞ ¼ ~ Aðr; zÞ:d~l. If the mutual inductance Lij between two coil sections is going to be calculated, then the average flux linking section i need to be calculated. This average flux is given by
R Ri R Zi /av g ¼
Ri
/ðr; zÞ dzdr i Ri Þ H i ðR Zi
ð3Þ
Knowing the average flux, the mutual inductance can be calculated using the following expression.
Lij ¼
Ni Nj /av g Ij
ð4Þ
White’s final expression for the mutual inductance between two coil segments is [13]:
Lij ¼ Lij0 þ 2Ni Nj ð1 mr Þl0 Rc
Z
1
0
I0 ðxRc Þ I1 ðxRc ÞFðxÞ dx mr þ ð1 mr ÞxRc I1 ðxRc ÞK 0 ðxRc Þ
ð5Þ
Fig. 3. Dimensions of two coils round an infinite core leg.
where
"
# Z xRi 1 2 FðxÞ ¼ xK 1 ðxÞdx x xðR i Ri Þ xRi xH i # " Z xRj xH i 1 2 sin xK 1 ðxÞdx 2 xHj xðR j Rj Þ xRj xH j sin cosðxdij Þ 2 1
ð6Þ
and where Lijo is the air-core inductance; mr = 1/lr is the relative reluctivity; and I0, I1, K0 and K1 are modified Bessel functions of first and second kind. Other parameters used in Eq. (6) are shown in Fig. 3. 3.2. Calculation of air-core inductances With respect to inductance definition, mutual inductance of two winding could be calculated from following equation.
Lij0 Ii Ij ¼
Z
~ Ai ~ J j dv
ð7Þ
~ Ai is the magnetic vector potential caused from winding i that in cylindrical coordinate system it is calculated as follows:
Ai ¼
l0 4p
Z
Z i
Zi
Z
R i
Ri
Z 2p 0
0
0
a/0 r0 dr d/0 dz J 1~ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 þ r 02 2rr 0 cosð/ /0 Þ þ ðz z0 Þ2
ð8Þ
In above equation, ðr; /; zÞ is the point in which the vector potential is calculated. The prime coordinates denote the source space (winding i). The values of integral limits are shown in Fig. 3. Vector ~ Jj a/ . After substituting values of ~ could be written as Jj~ Ai and ~ Jj in Eq. (7) and doing some simplifications (integrals by height are taken analytically), the final relation for the air-core mutual inductance could be written in a form of a triple integral as follows:
2568
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
Fig. 4. Comparing the inductances calculated by FEM and analytical method (a) air-core and (b) core leg.
Lij0 ¼
l0
Z
R i
Z
R j
Hi Hj Bwi Bwj Ri Rj Z p FðZi Zj Þ FðZ i Z j Þ FðZ i Z j Þ þ FðZ i Z j Þ dR1 dR2 d/ 0
ð9Þ where
FðZÞ ¼
R21 R22
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R21 þ R22 2R1 R2 cosð/Þ þ Z 2 R21 þ R22 2R1 R2 cosð/Þ
2
sin ð/Þ
ð10Þ
The integrals of the Eq. (9) are solved numerically. The integral by the angle is solved using Gauss–Legendre quadrature formula and the integrals by the elements radial dimensions are solved with the simpson’s 1/3 rule formula. 3.3. Comparing the calculated inductances with the FEM results In order to calculate the inductances by finite elements method, first all parts (HV disks/LV layers) are modeled in 2D-axisymmetric in a FEM tool. Then each part is defined as a stranded coil with specific turns. For calculation of the self and mutual inductances of each coil, current of I is applied to the desired coil while the current of other coils is defined as zero. After solution, the linkage flux (k) of each coil is determined. Because the system is linear, inductances could be calculated using the equation of L = k/I.
The values of calculated inductances for middle disk (disk 7) of HV winding by two methods in two cases of air-core and core leg (with core leg) are shown in Fig. 4. According to Fig. 4a, in the case of air-core, the results of two methods are completely close to each other. In the case of core leg since in the analytical method the height of core leg is assumed infinitely long, so for comparing results of two methods the height of core leg should be modeled long enough in the FE domain. As it is shown in Fig. 4b, by increasing the height of the modeled core leg, the FEM results becomes closer to the analytical method results. Inductances of HV winding disks calculated by Eq. (5) are shown in Fig. 5. 4. Capacitance calculation The most accurate method for calculating capacitance is using FEM. In order to calculate the capacitances between transformer windings using FEM, first all foils and all different dielectric materials are modeled in a FEM tool and for each set of different voltages applied to the conductors the total resultant energy is computed. Then the energy equations are solved simultaneously and the capacitances are obtained. By this method all the foils capacitances are calculated. But since all the capacitances are not actually used in the detailed model (because of negligible value) and also because of long solution time of this method it is better to apply substructuration principles to reduce the model order. For example with respect to construction of LV winding in dry transformers, the equation of plane capacitor could be used for calculation of lumped capacitance between two adjacent foils while the capacitance between two non-adjacent foils could be ignored because it is insignificant. The arrangement of conductors and dielectric materials such as air and resin and also the capacitance network is shown in Fig. 6. In the following, first the analytical method for calculation of the capacitances is described and then the results are compared with the corresponding numerical values obtained by FEM. 4.1. Calculation of capacitance between core and first layer of LV winding
Fig. 5. Inductances of HV winding disks calculated using analytical method.
With respect to large height of LV winding and ignoring the fringing effect of field in edges, we could calculate this capacity by means of the equation of cylinder capacitor. Since the core is grounded, this capacitance is actually the capacity of first layer of LV winding with the earth. Regarding that the 1st layer will create an equipotential surface in front of the core leg, other layers will
2569
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
LV foils
n n −1 n − 2
HV disks
i
4
3
2
1
+V CT
0 Volt Core leg
(n − 1) CD (i − 1)
Z CD
CD
0 Volt
CT
−V
Air Resin Fig. 6. The arrangement of conductors and dielectric materials between HV and LV windings.
have ignorable capacity with earth and we will ignore them in the detailed model. 4.2. Calculation of capacitance between LV winding layers With respect to the large height and internal diameter of LV winding and also the close distance between the LV winding layers, we could use the equation of plane capacitor to calculate the capacitance between two adjacent LV layers. 4.3. Calculation of series capacitance All of the disks in HV winding are wounded in one direction that is from inside to outside. By applying a voltage difference across a pair of disks terminals, the equipotential lines in the space between disks will be as shown in Fig. 7. As shown in Fig. 7, the equipotential lines are in the form of oblique lines.The usual method for calculating disk series capacitance is using the energy calculation method. If the energy of a disk, WT could be estimated then the disk series capacitance could be calculated from equation WT = 1/2CV2. In order to calculate the energy stored in a disk it is required to know the voltage distribution inside disk. In this approach the voltage is assumed to be evenly distributed within the disk. Assume that a voltage difference equal to 2 V is applied
Fig. 8. A pair of disks under applying a voltage difference equal to 2 V.
to terminals of a pair of disks (see Fig. 8). So the voltage across each disk will be equal to V. Calculation of the series capacitance of each disk is explained in the following. If the term CT denotes capacitance between touching turns and CD denotes capacitance between a turn of one disk and the corresponding turn of the other disk, CTR, the resultant inter-turn capacitance, is given by
C TR ¼
CT ðn 1Þ
ð11Þ
According to Fig. 8, the capacitance between conductor i and the zero-value equipotential surface and also the voltage across this capacitance are given by
C Di ¼
ðn 1Þ CD ði 1Þ
ð12Þ
V Di ¼
ði 1Þ V ðn 1Þ
ð13Þ
For the resultant inter-disk capacitance, CDR, following expression can be written. n 1X 1 C Di V 2Di ¼ C DR V 2 2 i¼2 2
ð14Þ
And so
C DR ¼
n CD 2
ð15Þ
The resultant series capacitance of each disk is given as the addition of the CTR and CDR. For disks which are adjacent to other disks from both down and up sides another term of CDR should be added to Cs. 4.4. Calculation of capacitances between HV and LV windings
Fig. 7. Equipotential lines between a pair of disks.
Since the last layer of LV will create an equipotential surface in front of other layers of LV, only the capacitance of last layer of LV to HV disks will be considered and the capacitance of other layers of LV to HV disks will be ignored, because they are very small. Therefore only the capacitances between the last layer of LV and the first layer of HV disks are calculated. Although calculation of capacitances between LV and HV windings using analytical method has
2570
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
not enough accuracy but it is possible to correct the analytical results by a correction factor obtained from experiments or numerical method solutions. For estimation of the capacitance between each disk of the HV and the last turn of the LV the following relation could be used.
Cw ¼ a
e 0 D m pH tresin =eresin þ t air =eair
ð16Þ
where Dm is average of the LV outer diameter and the HV inner diameter. tresin and tair are thicknesses of resin and air space respectively. H is height of HV disk and a is a correction factor which is different for each disk. It is found after FEM simulations of different dry-type transformers that the value of a is in the range of 1.1–1.4. From various FEM results following values are proposed for a: For two first upper disks and lower disk, a 1.3. For tap disks, a 1.35. For other disks, a 1.2.
tors are not meshed in the FEM method. Although this will reduce the number of elements but thickness of insulators in HV disks is too small and it is difficult to have very small elements inside each insulator. To overcome this problem, each two or more insulators (depending on required precision) are modeled as one insulator. Since in times less than 1 ls, the winding inductances have no effect on the initial voltage distribution (current in an inductance cannot be established instantaneously), the voltage distribution is predominantly decided by the capacitances. The initial voltage distributions based on capacitances calculated by two different methods of analytical and FE are shown in Fig. 9. The presented results in Fig. 9 are for a situation that the LV winding from both sides and the HV winding from its end are grounded. Comparing two curves of Fig. 9 it can be concluded that the proposed analytical method has enough accuracy for calculation of winding capacitances. 5. Losses calculation
4.5. Comparing the calculated capacitances with the FEM results For calculation of capacitances using FEM, the conductors which the capacitances related to them are not needed are defined as float conductors and then the capacitances between desired conductors are calculated by the method explained before (energy method). Hence, in addition to considering the effects of all conductors, there is no need to calculate all capacitances and so the solution time is decreased considerably. Since electric field in a conductor is zero and energy is stored only in dielectrics, conduc-
If the magnitude and direction of magnetic field around a conductor were determined, by knowing the dimensions of conductor, losses in the conductor could be calculated assuming that the amounts of Hy and Hx around the conductor are nearly constant. In fact, by using the Maxwell equations magnetic field inside conductor could be found and after calculation of current density by J ¼ r H, losses of conductor could be determined. Here, losses created in a conductor with the dimensions of Bo and Ho and the length of Lo (Fig. 10a) is calculated. Suppose that the field is only in the y axis direction and it changes along the x axis. It means,
! ! H ¼ HðxÞ j
ð17Þ
The diffusion equation for H can be written as 2
d H 2
dx
jxlrH ¼ 0
ð18Þ
By solving the Eq. (18), H inside the conductor will be as follows:
H ¼ Aemx þ Bemx
ð19Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffi In which m ¼ jxlr. By replacing the boundary conditions in Eq. (19), A and B are given by Bo
A¼ Fig. 9. Initial voltage distribution.
Bo
H1 em 2 H2 em 2 emBo emBo
~ ¼ HðxÞ~ ~ ¼ Hx ðyÞ~ Fig. 10. Estimation of eddy current losses in a conductor in two cases for leakage field. (a) H j and (b) H i þ Hy ðxÞ~ j.
ð20Þ
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574 Bo
B¼
2571
Bo
H1 em 2 H2 em 2 emBo emBo
ð21Þ
After determining H, J could be calculated as follows:
! @H ! ~ J¼ k ¼ mðAemx Bemx Þ k @x
ð22Þ
The losses produced in the conductor are equal to:
P¼
Z
RðdiÞ2 ¼
Lo
Z
r
H0 2
H0 2
Z
B0 2 B
J 2 dxdy
ð23Þ
20
In which r is conductivity of conductor. Now, suppose that the field around the conductor is in two x and y directions and its equation is as follows (Fig. 10b):
! ! H ¼ Hx ðyÞ~i þ Hy ðxÞ j
ð24Þ
In this case, the amount of H in each of two directions of x and y inside the conductor is calculated from: mx
mx
ð25Þ
Hx ¼ Cemy þ Demy
ð26Þ
Hy ¼ Ae
þ Be
The coefficients of A, B, C and D are calculated by using the boundary conditions. In this case, J is calculated from the following equation:
@Hy @Hx ! ~ k J¼ @x @y
ð27Þ
When the J is calculated, losses could be determined by Eq. (23). For using the above equations for calculation of winding losses it is required to know the magnetic field distribution inside winding. Consider a conductor carrying current It and with dimensions of ‘a’ and ‘b’ as shown in Fig. 11. According to Biot Savart equation, magnetic field intensity due to this conductor is given by [21],
y
P(x,y)
b
x
a Fig. 11. A conductor dimensions for calculation of magnetic field intensity components.
Fig. 12. Magnetic field in the z direction (axial field).
Fig. 13. Magnetic field in the r direction (radial field).
" 2 2 x þ 2a þ y þ 2b It 1 b Hy ðx; yÞ ¼ yþ ln 2 2 2 2pab 2 x 2a þ y þ 2b 2 ðx þ 2a Þ2 þ y 2b 1 b ln y 2 2 2 2 x 2a þ y 2b ! !! y þ 2b y 2b a
arctan arctan þ xþ 2 x þ 2a x þ 2a ! !!# y þ 2b y 2b a
arctan arctan x 2 x 2a x 2a
ð28Þ
For calculation of magnetic field in the x direction, a similar equation is used just with substituting terms of ‘x’ with ‘y’ and ‘a’ with ‘b’ in the above equation. Total magnetic field distribution inside a winding is obtained by adding the effect of all conductors. The magnetic field of dry transformer in directions of z and r is shown in Fig. 12 and Fig. 13 respectively. The field distribution in all disks except the end disks and also tap disks is approximately analogous. In the first and the last disks, radial component of the field is greater which leads to increase in losses. There is a similar condition for LV layers which the effect of field’s radial component is considerable. 5.1. Comparing the calculated resistances with the FEM results Since there are a lot of foils in each disk, modeling of all winding’s foils simultaneously in the FEM domain will be time consuming due to increasing in the number of elements. To overcome this problem, losses for each disk are calculated separately. To do so, only the desired disk is modeled as a set of foil conductors and the other disks are modeled as a solid non-conductive region with an equivalent ampere-turn density. Then a current equal to the winding current is applied to each foil of the desired disk. After solution, foils losses are obtained. The total losses of the disk are equal to sum of the losses of foils. Using the analytical formula makes it possible to estimate losses in high frequencies where the FEM solution will be time consuming due to increasing in the number of elements. If the eddy current losses are being analyzed by using FEM analysis, there should be at least one element in one skin-depth. Therefore the number of elements will increase and the solution especially in high frequencies will be time consuming. Nevertheless using analytical formula especially in frequencies above 200 kHZ lead to erroneous results. This is the frequency where the skin-depth is equal to the dimensions of conductors. The calculated resistances according to frequency using analytical formula and FEM for a disk of HV winding are shown in Fig. 14. As seen from FEM results the
2572
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
known input voltage. By multiplying X(jx) by H(jx) we can find the desired node voltage in frequency domain U(jx):
UðjxÞ ¼ XðjxÞ HðjxÞ
ð34Þ
As a last step, node voltages in time domain can be calculated by applying IFFT (Inverse Fast Fourier Transformation) as below:
uðtÞ ¼ IFFT½UðjxÞ
ð35Þ
8. Experiment results
Fig. 14. Calculated resistances for LV layers by analytical method, FEM and modified analytical method.
characteristic ‘‘knee’’ can be observed just below 200 kHZ where the analytical results are not good estimated above this frequency. In spite of this inaccuracy, little differences in output voltage waveforms were observed which could be ignored with a good approximation considering the reasonable short calculation time. But it is also possible to use FEM results to modification of analytical results for increasing the precision. It can be seen from Fig. 14 that the slope beneath the ‘‘knee’’ is proportional to the square of the frequency while the slope above the ‘‘knee’’ is proportional to the square-root of the frequency [20]. So in order to correct the analytical results at frequencies above the ‘‘knee’’, following equation could be used.
Rac ðf Þ ¼ 10ðlogðRac
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
FEM ð2e5ÞÞþ
logðf =2e5ÞÞ
;
f > 2e5 HZ
In order to find the impulse voltage distribution in transformers, some tests were performed on a high-voltage winding. This winding was manufactured in Iran-Transfo Company in such a way that there are tributaries from its every disk. Accessibility of outlets from each disk facilitated voltage measuring in every disk and so the impulse voltage distribution inside the winding could be found. In order to perform the tests, the high-voltage winding with a low-voltage winding were installed in the form of complete phase in the middle leg of a three-leg core. The voltage ratio in windings was 20/0.4 kV. The test object is shown in Fig. 15. The impulse and FRA tests were done on the test object. The impulse tests include full and chopped-wave lightning tests. The measuring circuit of impulse tests is shown in Fig. 16. In the full lightning test, the impulse voltage of 1.2/50 ls was applied on
ð29Þ
The calculated resistances using modified analytical method are shown in Fig. 14 which are in good agreement with the FEM results. 6. System description The system description can be used to investigate internal resonances of the transformer winding, or to determine transfer-functions between terminals. The nodal admittance matrix (without considering the inductive branches) and the branch impedance matrix are as follows:
Y N ðjxÞ ¼ GðxÞ þ jxCðxÞ
ð30Þ
Z B ðjxÞ ¼ RðxÞ þ jxLðxÞ
ð31Þ
Fig. 15. Test object.
The branch-matrix ZB(jx) is then transformed to a nodal form by means of a transformation-matrix:
Y B ðjxÞ ¼ AZ B ðjxÞ1 AT
Y SYS ðjxÞ ¼ Y N ðjxÞ þ Y B ðjxÞ
ð33Þ
7. Calculation of node voltages At the first step, H(jx) (the voltage transfer function of each node to input) is calculated using detailed model. At the second step, the input impulse voltage in frequency domain X(jx) can be obtained by applying the FFT (Fast Fourier Transformation) to the
Disk 1
Disk 2
Disk 2
Disk 3
Disk 3
Disk 4
Disk 4
Disk 5
ð32Þ
The transformation-matrix A, describes the relation between nodal currents and branch currents. The system admittance matrix is then given as:
Disk 1
LV
Disk 6
LV
Disk 7 Disk 8
HV
Disk 5 Disk 6 Disk 7 Disk 8
Disk 9
Disk 9
Disk 10
Disk 10
Disk 11
Disk 11
Disk 12
Disk 12
Disk 13
Disk 13
(a)
HV
(b)
Fig. 16. Measuring circuit (a) impulse test and (b) FRA test.
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
2573
Terminal conditions in all simulations are similar to the test arrangement FRA analysis results for disk 3, disk 4 and disk 5 are shown in Figs. 19–21 and compare with test results. The impulse voltage results from test and simulation for full and chopped lightning impulses are shown in Fig. 22 and Fig. 23 respectively. Comparison between simulation results and experimental results shows that the proposed method is suitable for high frequency modeling of dry-type transformers.
Fig. 17. Measured disk voltages to ground at full lightning impulse.
Fig. 19. FRA results for disk 3.
Fig. 18. Measured disk voltages to ground at chopped lightning impulse.
HV winding input and then the output voltage of disks were measured and stored by means of a digital oscilloscope. These tests for chopped lightning impulse also repeated. Measurement results for full and chopped impulse tests are shown in Fig. 17 and Fig. 18 respectively. In FRA test, a sinusoidal voltage with various frequencies is applied on transformer input and in each frequency, the amplitude and phase of desired output sinusoidal voltage were measured and then by using the amplitude and phases of input and output voltages, the amplitude and phase of transfer function in that frequency could be calculated. In order to perform FRA test, a FRA analyzer device which was fabricated in OMICRON Company was used. This device has the capability of sweeping the input sinusoidal voltage from 1 kHZ to 5 MHZ and also it can measure the desired output voltage. This device has a connection with computer and specific software can instantly calculate and draw output to input transfer function. The input signal was applied to HV winding input and then the frequency response of each disk was measured. In the FRA tests, as shown in Fig. 14, output signals were measured from the top of a 50 ohm shunt resistance which was grounded.
Fig. 20. FRA results for disk 4.
9. Simulation results To verification the proposed method, parameters of detailed model were calculated for the test object. After establishing the model it was employed in MATLAB and FRA and impulse voltage analysis were implemented.
Fig. 21. FRA results for disk 5.
2574
M. Eslamian et al. / Energy Conversion and Management 52 (2011) 2565–2574
Fig. 22. Simulated and measured disk voltages to ground at full lightning impulse.
Fig. 23. Simulated and measured disk voltages to ground at chopped lightning impulse.
10. Conclusion In this paper high frequency behavior of cast resin dry-type transformers was simulated. First a detailed model proposed for these transformers and then parameters of the model were calculated using analytical formula and then compared with the corresponding FEM results. In order to verify the model’s accuracy a lab transformer was constructed and FRA and impulse tests were implemented. Comparison between simulation and experimental results shows that the proposed method is well suited for modeling the transient behavior of cast resin dry-type transformers. References [1] Lewis TJ. The transient behavior of ladder networks of the type representing transformer and machine windings. IEE Proc 1954;101(pt. II):541–53. [2] McWhirter JH, Fahrnkopf CD, Steele JH. Determination of impulse stresses within transformer windings by computers. AIEE Trans 1957;75(pt. III):1267–79. [3] Dent BM, Hartill ER, Miles JG. A method of analysis of transformer impulse voltage distribution using a digital computer. IEE Proc 1958;105(pt. A):445–59. [4] Okuyama KA. A numerical analysis of impulse voltage distribution in transformer windings. Electr Eng Jpn 1967;87(no. 1):80–8.
[5] Stein M. A study of the initial surge distribution in concentric transformer windings. AIEE Trans 1964:877–92. [6] Kawaguchi Y. Calculation of the circuit constants for computing internal oscillating voltage in transformer windings. Electr Eng Jpn 1969;89(no. 3):44–53. [7] Miki A, Hosoya T, Okuyama K. A calculation method for impulse voltage distribution and transferred voltage in transformer windings. IEEE Trans Power App Syst 1978;PAS-97:930–9. [8] De Leon F, Semlyen A. Efficient calculation of elementary parameters of transformers. IEEE Trans Power Deliv 1992;7(no. 1):376–83. [9] De Leon F, Semlyen A. Complete transformer model for electromagnetic transients. IEEE Trans Power Deliv 1994;9(no. 1):231–9. [10] Rabins L. Transformer reactance calculations with digital computers. AIEE Trans 1956;75:261–7. [11] Fergestad PI, Henriksen T. Inductances for the calculation of transient oscillations in transformers. IEEE Trans Power App Syst 1974;PAS-93:510–7. [12] Fergestad PI, Henriksen T. Transient oscillations in multi-winding transformers. IEEE Trans Power App Syst 1974;PAS-93:500–9. [13] James H. Harlow, electric power transformer engineering, transient-voltage response (chapter 20). Robert C. Degeneff, CRC Press; 2004. [14] Wilcox DJ, Conlon M, Hurley WG. Calculation of self and mutual impedances for coils on ferromagnetic cores. IEE Proc Sci Meas Technol 1988;135(pt. A, no. 7):470–6. [15] Wilcox DJ, Hurley WG, Conlon M. Calculation of self and mutual impedances between sections of transformer windings. IEE Proc Gener Trans Distrib 1989;136(pt. C, no. 5):308–14. [16] Yamashita H, Cingoski V, Nakamae E, Namera A, Kitamura H. Design improvements on graded insulation of power transformers using transient electric field analysis and visualization technique. IEEE Trans Energy Convers 1999;14(no. 4):1379–84. [17] Miri AM, Riegel NA, Kuhner A. Finite element models for the computation of the transient potential and field distribution in the winding system of high voltage power transformers. IEE Conf High Voltage Eng 1999;2(467):39–42. [18] Azzouz Z, Foggia A, Pierrat L, Meunier G. 3D finite element computation of the high frequency parameters of power transformer windings. IEEE Trans Magn 1993;29(no. 2):1407–10. [19] Moreau O, Michel R, Chevalier T, Meunier G, Joan M, Delcroix JB. 3-D high frequency computation of transformer R,L parameters. IEEE Trans Magn 2005;41(no. 5):1364–7. [20] Bjerkan E, Høidalen HK. High frequency FEM-based power transformer modeling: investigation of internal stresses due to network-initiated overvoltages. Electr Power Syst Res 2007;77(11):1483–9. [21] Binns KJ, Lawrenson PJ, Trowbridge CW. The analytical and numerical solution of electric and magnetic fields. John Wiley and Sons; 1994. Morteza Eslamian was born in Zanjan, Iran, in November 1984. He received the B.S. degree in electrical engineering in 2006 from Zanjan University and the M.Sc. degree in electrical engineering in 2009 from Amirkabir University of Technology, Tehran, Iran, where he is currently pursuing the Ph.D. degree. His research interests include the applications of neural network and genetic-algorithm (GA) in power systems analysis and finite element method (FEM) in transformer high frequency and thermal modeling. He served in Iran-Transfo Co. high voltage lab from 2005 to 2006. He is a member of Iran-Transfo Co. R&D from 2005. Behrooz Vahidi was born in Abadan, Iran, in 1953. He received the B.S. degree in electrical engineering from Sharif University of Technology, Tehran, Iran, in 1980, the M.S. degree in electrical engineering from Amirkabir University of Technology, Tehran, Iran, in 1989, and the Ph.D. degree in electrical engineering from the University of Manchester Institute of Science and Technology, Manchester, UK, in 1997. From 1980 to 1986, he worked in the field of high voltage in industry as Chief Engineer. From 1989 to the present, he has been with the Department of Electrical Engineering of Amirkabir University of Technology, where he is now a Professor. He is IEEE senior member. His main fields of research are high voltage, electrical insulation, power system transient, lightning protection, and pulse power technology. He has authored and coauthored more than 200 papers and five books on high-voltage engineering and power system. Seyed Hossein Hosseinian was born in 1961 in Iran. He received both the B.Sc. and M.Sc. degrees from the Electrical Eng. Dept. of Amirkabir University of technology, Iran, in 1985, and 1988, respectively, and PhD degree in Electrical Engineering Dept, university of Newcastle England, 1995. At the present, he is the Associate Professor of Electrical engineering Department in Amirkabir University of technology (AUT). His especial fields of interest include transient in power systems, Power Quality, Restructuring and Deregulation in power systems. He is the author of four books in the field of power systems. He is also the author and the coauthor of over one hundred technical papers.