Unification of anisotropy and FEM-BE models for distribution transformer optimization

JOURNAL OF OPTOELECTRONICS AND ADVANCED MATERIALS Vol. 10, No. 5, May 2008, p. 1143 - 1148 Unification of anisotropy and FEM-BE models for distribution transformer optimization THEMISTOKLIS KEFALAS*, MARINA TSILI, ANTONIOS KLADAS Laboratory of Electrical Machines, Electric Power Division, Faculty of Electrical and Computer Engineering, National Technical University of Athens, 9 Heroon Polytechniou Street, Athens 15780, Greece The paper presents a three-dimensional finite element (3D FEM) anisotropy model, based on a particular scalar potential formulation, for the no load loss evaluation of wound core shell type distribution transformers. The specific 3D FEM anisotropy model is combined with a hybrid finite element-boundary element (FEM-BE) model, used for the calculation of the transformer’s short circuit impedance, and various optimization algorithms in order to minimize the total owing cost (TOC) of a distribution transformer. (Received March 13, 2008; accepted May 5, 2008) Keywords: Transformer, Finite element anisotropy model, Hybrid FEM-BE model 1. Introduction The design optimization of a wound core shell type distribution transformer consists in minimizing the transformer's total owing cost (TOC), where TOC is defined as the first cost of the transformer plus the calculated present value (PV) of its future losses [1]. The transformer manufacturer must also take into account the transformer's ratings, design constraints, and specifications imposed by the electric power supplier [2]. Key specification parameters that must be met by the manufacturer are the no load loss and the short circuit impedance. The evaluation of the aforementioned parameters is based most of the times on analytical methods [2]. However a simple analytical method for the computation of the short-circuit impedance and the no load loss does not suffice especially in the case of transformers that do not fit the standardized large-scale constructions [3]. As a result the manufacturer is forced to adopt a safety margin for the no load losses [4], [5], and to install magnetic shunts in order to correct within the desired limits the short circuit impedance [6]. In turn those actions result in the increase of the manufacturing cost and losses. The paper addresses the latter problems by introducing an appropriate three-dimensional finite element (3D FEM) anisotropy model used for the accurate prediction of the transformer's no load losses. The specific 3D FEM anisotropy model is combined with a hybrid finite element-boundary element (FEM-BE) model [3], [6], used for the evaluation of the short circuit impedance, through the process of minimization of the TOC. 2. 3D FEM formulation Solving 3D magnetostatic problems by the finite element (FE) method based on the magnetic vector potential (MVP) is quite laborious as there are three unknown components to be determined at each node of the FE mesh [7]. This is the most obvious drawback of the MVP formulation but it is not the only one. The current carrying regions (coils) must be discretized carefully in order to take their effect into account accurately, and the solution has be found to be erroneous when the normal component of the MVP is significant at the interface between regions of different permeability, e.g. at the interface between iron and air [7], [8]. By using edge elements instead of nodal, the continuity of the tangential component of the MVP is satisfied without constraining the normal components on element boundaries. However the resulting formulation is more complicated and edge elements are numerically less stable than nodal elements [7]. On the other hand magnetic scalar potential (MSP) based FE formulations are advantageous in terms of computational effort, as there is only one degree of freedom to be evaluated for every node of the FE mesh. Despite this obvious advantage, many problems arise with the use of the MSP like inherent difficulties in handling 3D volume current distributions, cancellation errors in the iron domain, difficulties in modeling multiply connected iron regions, and multiply valued potentials [7]. A number of scalar potential formulations have been developed to address some of the latter problems like the difference potential (DP) and the total scalar potential (TP) formulation. Those early formulations have been united and extended by the general potential (GP) formulation [8]. According to the GP formulation the magnetic field intensity H is sought in the following form H = H g + ∇Φ g (1) 1144 Themistoklis Kefalas, Marina Tsili, Antonios Kladas where H g is an initial guess magnetic field and Φg is the general potential. If H g satisfies Ampere's law and its absolute value is much larger than that of ∇Φ g then the ∇ ⋅ (μ ⋅ (H g + ∇Φ g )) = 0. solution of the problem can be found according to (2) What remains is the evaluation of a suitable guess magnetic field H g with the following three-step scheme proposed in [8]. 1. First compute the guess magnetic field H gi , in conductors shape with little computational effort. The problem's solution is obtained by discretizing (6) that ensures the total's field solenoidality ∇ ⋅ (μ ⋅ (K − ∇Φ )) = 0. (6) The previous formulation does not suffer from cancellation errors, it satisfies Ampere's law, and consequently it is applicable to multiply connected iron regions. Simplicity and computational efficiency are its main advantages in contrast with the GP formulation, rendering the specific MSP formulation ideal for providing the solution to a stochastic optimization algorithm. the iron domain by satisfying the following conditions. ∇ ⋅ (μ ⋅ H gi ) = 0, n ⋅ (μ ⋅ H gi ) = 0. 2. 3. Anisotropy FEM model (3) Then compute the guess magnetic field H go , in the air domain by satisfying the following conditions. ∇ ⋅ (μ ⋅ H go ) = 0, ∇ × H go = J o , n × (H gi − H go ) = 0. (4) 3. Calculate the general potential Φg over the whole domain by satisfying (2). Even though the GP formulation covers successfully most of the practical engineering problems there are two disadvantages that might make the particular formulation not suitable for using it in conjunction with optimization algorithms, especially stochastic. The first one is the complexity of the GP formulation since a three-step procedure is required to obtain the solution. The second disadvantage is that for a nonlinear analysis several iterations are required during the execution not only of the third step but also of the first step for the calculation of the scalar potential in the iron region, Φgi . These disadvantages lead to increasing computational effort and time. In this paper, it is proposed to use a particular scalar potential formulation necessitating no prior source field calculation, presented in [7]. According to the specific formulation H is partitioned as follows H = K − ∇Φ By considering the iron-laminated material as homogeneous and anisotropic media at the level of finite elements an accurate representation of the core material is achieved [9]. An elliptic anisotropy model is best suited for the wound core transformer in contrast with the stack core transformer [10]. The specific model is based on the assumption that the field intensity H has an elliptic trajectory for the modulus of the flux density constant. Therefore, if μ p is the magnetic permeability tangential to the electrical steel rolling direction, μ q is the permeability normal to the lamination, and r is the ratio of the ellipse semi-axes then μ q = rμ p , 0 < r < 1 . (7) The permeability tensor μ in the global coordinate system is given by (8), where R is the rotation matrix, and μ F is the permeability tensor in the local coordinate system. μ = R −1 μ F R (8) The accurate evaluation of the flux density distribution with the specific FE method is used in conjunction with the experimentally determined specific core losses, for the evaluation of the wound core no load loss. 4. FEM-BE model (5) where Φ is a scalar potential extended all over the solution domain and K is a fictitious field distribution that satisfies the following three conditions. 1. K is confined in a simply connected subdomain comprising the coil region. 2. ∇ × K = J in the coil region and ∇ × K = 0 outside the coil domain. 3. Lastly K is perpendicular on the boundary of the subdomain comprising the coil region. The above conditions constitute a minimum enabling to simulate the coil domain. The distribution of K is easily determined analytically or numerically by the During the short-circuit test, the electromagnetic field is not confined only in the wound cores and the conductors but it expands over extensive parts of air between the transformer’s active part and tank walls. Therefore for the evaluation of the short-circuit impedance expect from the accurate representation of the low voltage (LV) and high voltage (HV) winding geometry, a dense discretization of the air surrounding the transformer’s active part is needed. This in turn leads in considerable computation cost in the case of the FE method. We have addressed the specific problem by developing a hybrid FEM-BE technique [3], [6]. The BE method uses the integral form of magnetic Unification of anisotropy and FEM-BE models for distribution transformer optimization 1145 field equations, that is mathematically equivalent to the original partial differential equation. The boundary integral equation corresponding to Laplace equation is of the form ∂Φ(s ′)⎤ ∂G (s ′,s ) ⎡ ds ′ = 0 c(s ) Φ(s ) + Φ ⎢Φ(s ) − G (s ′,s ) Γ ∂n ′ ⎥⎦ ∂n ⎣ (9) where s is the observation point, s ′ is the boundary Γ coordinate, n′ is the unit normal, and G is the fundamental solution of Laplace equation in free space. The BE method discretizes only the boundaries of the domain, in contrast to the FE method which discretizes the whole domain. This makes the BE method suitable for geometries with extensive parts of air. The use of the hybrid FEM-BE technique results in the accurate evaluation of the magnetic field in the transformer’s active part and in the air region surrounding it, with low computational effort comparing with the FE method. Fig. 1 depicts the 3D one-phase model of the three-phase wound core distribution transformer, consisting of the LV and HV windings and the iron cores. The hybrid FEM-BE model is divided in two subdomains, the transformer’s active part, FEM subdomain, represented by a tetrahedral FE mesh, and the subdomain between the active part and the transformer’s tank walls, BE subdomain, represented by a triangular mesh of its boundaries. Fig. 1. Geometry of the FEM-BE model. 5. Experimental verification of the 3D FEM models Fig. 2. Experimental setup. Fig. 3. Detail of the wound core, the excitation coil, and the search coils. Two turn search coils wound around the core’s sheets, were employed for determining the flux distribution along the core’s limb and corner for different magnetization levels. The voltages induced in the search coils were measured by connecting their terminals directly to the DAQ card’s voltage differential inputs. Fig. 3 illustrates a detail of the wound core, the excitation coil, and the search coils. The manipulation of the data captured by the DAQ card was performed by virtual instruments (VI) created with the use of LabView. Fig. 4 illustrates a perspective view of the geometry, and the flux density vector plot of a wound core, obtained by the 3D FEM anisotropy model. The core is built of high permeability magnetic steel, and the mean flux density used for the FE analysis is 1.57 T. Figs. 5 and 6 display the flux density vector, and nodal plot respectively, of the upper front part of the same wound core, for a mean flux density of 1.66 T. The 3D FEM anisotropy model presented in Sections 2 and 3 was used for the numerical analysis of various one-phase, and three-phase wound core transformers. The computed no load loss value, and flux density distribution at the core’s limb and corner, were compared to the measured ones. Fig. 2 shows the experimental setup used for the no load loss and local field measurements of onephase wound cores. It consists of a PC, a National Instruments (NI) 6143 data acquisition (DAQ) card, an active differential voltage probe and a current probe based on the Hall Effect. For the magnetization of the one-phase wound core, a 23 turn excitation coil and a 230 V, 50 Hz, one-phase autotransformer were used. Fig. 4. Perspective view of the geometry and the flux density vector plot of a one-phase wound core (B = 1.57 T). 1146 Themistoklis Kefalas, Marina Tsili, Antonios Kladas Fig. 5. Vector plot of the wound core flux density distribution during no load test (B = 1.66 T). Fig. 8. Flux density magnitude distribution along line CD. Fig. 6. Nodal plot of the wound core flux density distribution during no load test (B = 1.66 T). The computed flux density distribution across the core’s limb, line AB of Fig. 4, and corner, line CD of Fig. 4, for two different magnetization levels, 1.57 T and 1.66 T, are compared with the measured ones in Figs. 7 to 10. The experimentally defined curves show that the flux density distribution along the limb is different from the flux density distribution along the corner. In both cases the flux density is low at the most inner steel sheets, and then it increases. However the gradient of the flux density is much steeper in the case of the flux density distribution across line AB. The computed results, shown in Figs. 7 to 10, verify that the flux density distribution along lines AB and CD is different. Nevertheless, due to the fact that the FE analysis does not take into consideration the stress induced by manufacturing processes, the error of the computed flux density is greater in the inner sheets than in the outer sheets, especially along line AB as can be seen from Figs. 7 and 9. Fig. 9. Flux density magnitude distribution along line AB. Fig. 10. Flux density magnitude distribution along line CD. The error between calculated and measured no load loss, for a number of one-phase wound cores, was most of the times underestimated and less than 5 %. However this is not the case for the three-phase wound core transformer where the evaluation of the no load loss calls for multiple magnetostatic analysis. Table 1 summarizes the computed and measured no load loss of three-phase transformers of different apparent power and working induction ratings. The error between computed and measured no load loss is underestimated and also overestimated in some cases, whereas the error exceeds in one case 5 %. Fig. 7. Flux density magnitude distribution along line AB. 1147 Unification of anisotropy and FEM-BE models for distribution transformer optimization Table 1. Comparison of the computed and measured no load loss values, of three-phase wound core distribution transformers of various apparent power and working induction ratings. Apparent power (kVA) 400 630 1,600 250 1,000 630 Mean flux density (Tesla) 1.35 1.37 1.55 1.66 1.72 1.75 Computed no load loss (W) Measured no load loss (W) Error (%) 553 790 1,642 446 1,277 1,099 568 830 1,620 430 1,350 1,130 2.64 4.82 -1.36 -3.72 5.41 2.74 The results obtained by the hybrid FEM-BE technique were validated for a number of three-phase distribution transformers through local leakage field measurements by a Hall Effect probe during the short-circuit test. Also the difference between the measured and the computed short circuit impedance is in most cases less than 2.7 %. Further details on the experimental verification of the hybrid FEM-BE technique may be sought in [3] and [6]. 6. Wound core transformer design optimization procedure The unification of the 3D FEM anisotropy and the FEM-BE models is achieved by the minimization of the TOC described by (10), where C Fe and CCu are the magnetic steel and the winding material unit cost ($/Kg), M Fe and M Cu are the mass of the magnetic steel and the winding material (kg), SM is the sales margin, and PNLL , PLL are the no load and load loss (W) respectively. The A factor and B factor ($/W) are the PV of 1 W of no load loss and load loss respectively, over the life of the transformer [1]. The minimization of (10) is subject to the three inequality constraints of (11) according to IEC 60076-1 spec [5], where U k is the short circuit impedance (%), PNLL , PLLspec , and U kspec are the specified by international technical specifications [5], values of no load loss, load loss, and short circuit impedance. It is also subject to the apparent power constraint, the induced voltage constraint, the temperature rise constraint and the no load current constraint [2], [5]. TOC = (C Fe M Fe + CCu M Cu ) / SM + A factor PNLL + B factor PLL (10) spec PNLL < 1.15 PNLL , PLL < 1.15 PLLspec , U k < 1.1U kspec HV winding, cs LV and cs HV respectively. During the optimization process, the evaluation of PNLL and U k is performed by the 3D FEM anisotropy model of Section 3, and the FEM-BE model of Section 4. The evaluation of M Fe , M Cu , and PLL is based on simple analytical relationships [2], [5]. Also the thermal calculation of the transformer is realized through the number of the winding’s cooling ducts, shown in Fig. 5, and the calculation of the cross-sectional area of the conductors is implemented from the current density [5]. The proposed optimization procedure was applied to a 160 kVA, 20 / 0.4 kV, three-phase distribution transformer, consisting from two outer wound cores with design parameters as illustrated in Fig. 4, and two inner cores, with a yoke length F 2 twice as large as that of the outer cores ( F 2 = 2 ⋅ F1 ). The current density chosen for the LV and HV winding is 3.2 A / mm2 and 3.7 A / mm2 respectively, and the specified operating parameters are spec PNLL = 425 W, PLLspec = 2,350 W, and U kspec = 4 %. The specific optimization problem presents multiple optima and it was tackled by stochastic optimization algorithms. In [2] the geometric programming was applied, whereas in [5] a heuristic solution was chosen. In the present paper two well known stochastic algorithms were tested the genetic algorithm (GA) and the simulated annealing (SA). Table 2 summarizes the transformer’s parameters and the solution obtained by the SA algorithm, i.e. the most effective of the algorithms tested. The GA was proven to be more computationally expensive as it required a greater number of objective function evaluations, and also the final TOC value was 0.51 % higher than that obtained by the SA algorithm. Table 2. Three-phase wound core distribution transformer parameters and solution output. Parameter Apparent power (kVA) LV coil voltage (V) HV coil voltage (V) Connection of LV coil Connection of HV coil Frequency (Hz) Value 160 Solution G (mm) Value 216 400 D (mm) 190 20.000 B (T) 1.607 Y Np 30 D csLV (mm2) 72.17 50 HV conductor diameter (mm) 1.04 (11) The independent variables chosen for the solution of the specific optimization problem are G , D , depicted in Fig. 4, the flux density B , the number of turns of the LV winding N p , and the cross-sectional area of the LV and 7. Conclusion A 3D FEM anisotropy model for the no load loss evaluation of one-phase and three-phase wound core transformers was developed, based on a computationally efficient scalar potential formulation. The accuracy of the 1148 Themistoklis Kefalas, Marina Tsili, Antonios Kladas specific model was validated by iron loss and local field measurements. It was also combined with a hybrid 3D FEM-BE model, for the transformer design optimization. The accuracy and the computational efficiency of the proposed models are of great importance for transformer manufacturers as it will enable them to design transformers that do not fit the standardized large-scale constructions without the need to adopt safety margins. Acknowledgment This work was supported in part by the General Secretariat for Research and Technology of Greece under PENED Grant No. 03ED045. References [1] S. Y. Merritt, S. D. Chaitkin, IEEE Industry Applications Magazine, 9(6), 21 (2003). [2] R. A. Jabr, IEEE Trans. Magnetics 41(11), 4261 (2005). [3] M. Tsili, A. Kladas, P. Georgilakis, A. Souflaris, D. Paparigas, Journal of Materials Processing Technology 161, 320 (2005). [4] P. Georgilakis, N. Hatziargyriou, D. Paparigas, IEEE Computer Applications in Power 12(4), 41 (1999). [5] P. S. Georgilakis, M. A. Tsili, A. T. Souflaris, “A heuristic solution to the transformer manufacturing cost optimization problem,” Journal of Materials Processing Technology, 181, 260 (2007). [6] M. A. Tsili, A. G. Kladas, P. S. Georgilakis, A. T. Souflaris, D. G. Paparigas, IEEE Trans. Magnetics 41(5), 1776 (2005). [7] A. Kladas, J. Tegopoulos, IEEE Trans. Magnetics 28, 1103 (1992). [8] M. Gyimesi, D. Lavers, T. Pawlak, D. Ostergaard, IEEE Trans. Magnetics, 29(2), 1345 (1993). [9] T. D. Kefalas, P. S. Georgilakis, A. G. Kladas, A. T. Souflaris, D. G. Paparigas, IEEE Trans. Magnetics, to be published. [10] H. V. Sande, T. Boonen, I. Podoleanu, F. Henrotte, K. Hameyer, IEEE Trans. Magnetics 40(2), 850 (2004). ____________________________ * Corresponding author: [email protected]