Leakage Inductances of Transformers at Arbitrarily Located Windings

: The article presents the calculation of the leakage inductance in power transformers. As a rule, the leakage ﬂux in the transformer window is represented by the short-circuit inductance, which a ﬀ ects the short-circuit voltage, and this is a very important factor for power transformers. This inductance reﬂects the typical windings of power transformers well, but is insu ﬃ cient for special transformers or in any case of the internal asymmetry of windings. This paper presents a methodology for calculations of the self- and mutual-leakage inductances for windings arbitrarily located in the air window. It is based on the 2D approach for analyzing the stray ﬁeld in the air zone only, using discrete partial di ﬀ erential operators. That methodology is veriﬁed with the ﬁnite element method tested on real transformer data.


Introduction
In the long history of research on power transformers, countless books have been published, for example [1][2][3][4][5][6][7]. Theoretically, basic construction and operational problems were solved there. In addition to single-phase and three-phase transformers, new types of special power transformers for power system control, power electronics, and electric and traction have been designed. Two approaches are used to analyze the properties of transformers. The first is based on Maxwell's field equations. The second is based on Kirchhoff's equations with mutual-inductances. The field approach is mainly used for transformer design problems, the second for operational problems when the transformer is part of the system. The field approach uses the design data of magnetic circuits and windings, the circuital approach requires relations between the linked fluxes of lumped inductances of the transformer windings and the winding currents. Due to the magnetic non-linearity of the transformer ferromagnetic core, these relations are non-linear, which can only be determined by solving Maxwell's equations [8][9][10][11][12][13][14][15]. Nowadays, software packages allow us to combine the field and the circuit approaches and calculate output values like currents or voltages at transient or steady-states that are interesting for engineers [16]. However, this is costly and time-consuming and seems inappropriate for solving operational problems in systems with power transformers. Therefore, simplified methods for calculating the lumped parameters of a transformer are often used [17][18][19][20][21][22].
To overcome the difficulties of 2D or 3D non-linear magnetic field computation in the whole transformer magnetic circuit, the relations between linked fluxes ψ n and currents i n of elementary The above two matrices can be calculated separately using a simplified method. The first matrix has inductances, due to fluxes in the core, and the second one results from the leakage field in the transformer windows. Inductances L µ n,m (i) are non-linear, taking into account the saturation of the ferromagnetic core, and may depend on all winding currents. This matrix can be determined using circuit representation of ferromagnetic core with non-linear lumped parameters, omitting fluxes in the air. The second matrix, containing the leakage inductances, can be found taking into account only the window zone and assuming that µ Fe → ∞. The leakage inductances L σ n,m are constant because the window zone is magnetically linear.
For symmetrically designed power transformers, with magnetic circuits and winding locations sketched in Figure 1, the matrices in (1) take extremely simple forms, which are well known from books and textbooks. Commonly, those matrices are determined by three lumped parameters: Non-linear magnetizing inductance L µ (i µ ) and two leakage inductances of primary L σ p and secondary winding L σ s . The inductance L µ (i µ ) represents relation between the linked flux in the limb and the magnetizing current. The sum of leakages inductances determines the short-circuit inductance L sc , and consequently, the short-circuit voltage. It is a very important technical parameter of power transformers. The short-circuit inductance is calculated from the magnetic field in the transformer window with two windings on the same limb, assuming that the ampere-turn compensates and the magnetic permeability of the ferromagnetic core is sufficiently high. The field analysis in the 1D transformer window model provides rather simple formulae for the short circuit inductance. This can be found in books on transformers, and is still recommended for designers. The leakage inductances L σ p and L σ s are not calculated, but estimated based on the short circuit inductance [19,20], often just by dividing the short-circuit inductance by 2. This approach is sufficient for symmetrically located windings, when they cover, more or less, the limbs, and the Rogowski factor can correct the length of the magnetic field path in the air. Figure 2a illustrates such a case. Figure 2b shows the position of windings when HV windings are modified, for voltage regulation purposes, by eliminating some numbers of turns in the middle part. Figure 3a presents a case when a short circuit in the HV winding in phase 'A' arise. For the cases 'b' and 'c' leakage inductances cannot be estimated from 1D field analysis.  μ  μ  μ  σ  σ  σ  1  1  1,1  1,2  1,N  1,1  1,2  1,N  μ  μ  σ  σ  2  2  2,2  2,N  2,2 The above two matrices can be calculated separately using a simplified method. The first matrix has inductances, due to fluxes in the core, and the second one results from the leakage field in the transformer windows. Inductances μ , L ( ) n m i are non-linear, taking into account the saturation of the ferromagnetic core, and may depend on all winding currents. This matrix can be determined using circuit representation of ferromagnetic core with non-linear lumped parameters, omitting fluxes in the air. The second matrix, containing the leakage inductances, can be found taking into account only the window zone and assuming that Fe μ → ∞ . The leakage inductances σ , L n m are constant because the window zone is magnetically linear.
For symmetrically designed power transformers, with magnetic circuits and winding locations sketched in Figure 1, the matrices in (1) take extremely simple forms, which are well known from books and textbooks. Commonly, those matrices are determined by three lumped parameters: parameter of power transformers. The short-circuit inductance is calculated from the magnetic field in the transformer window with two windings on the same limb, assuming that the ampere-turn compensates and the magnetic permeability of the ferromagnetic core is sufficiently high. The field analysis in the 1D transformer window model provides rather simple formulae for the short circuit inductance. This can be found in books on transformers, and is still recommended for designers. The leakage inductances σ p L and σ s L are not calculated, but estimated based on the short circuit inductance [19,20], often just by dividing the short-circuit inductance by 2. This approach is sufficient for symmetrically located windings, when they cover, more or less, the limbs, and the Rogowski factor can correct the length of the magnetic field path in the air. Figure 2a illustrates such a case. Figure 2b shows the position of windings when HV windings are modified, for voltage regulation purposes, by eliminating some numbers of turns in the middle part. Figure 3a presents a case when a short circuit in the HV winding in phase 'A' arise. For the cases 'b' and 'c' leakage inductances cannot be estimated from 1D field analysis.     Rectifier transformers for multi-phase voltages and currents [7] can have windings that are much shorter than the limb's height, and the Rogowski factor cannot be applied. Figure 2a shows an example of the location of the windings in the transformer window for Yyd 3/6-phase transformers for 12 pulses rectifiers and Figure 2b for Yzyd 3/12-phase transformers for 24 pulses rectifiers. The position of the windings loses the characteristic symmetry for classic three-phase power transformers, see Figure 1. At least the 2D field is necessary to find leakage inductances.
An analytical solution of the magnetic vector potential in transformer air windows with arbitrarily located balanced pairs of windings is presented in [21]. The zero Neumann boundary conditions can be stated for such cases on the ferromagnetic boundaries, which means that the magnetizing current is omitted. We can then calculate the short circuit inductance for any pair of windings. However, this approach cannot be used to determine the leakage inductances of individual windings and the mutual-inductances caused by the stray field. To find these inductances, 2D or 3D field calculations are typically needed, including the entire magnetic circuit with the ferromagnetic core and the air zone. [11][12][13][14][15]22]. This paper presents a 2D numerical approach for determining all leakage inductances σ , L n m , both self, as well as mutual, for arbitrarily located windings considering the magnetic field in the air window only, when surrounded by ferromagnetic walls with high permeability. Difficulties in    Rectifier transformers for multi-phase voltages and currents [7] can have windings that are much shorter than the limb's height, and the Rogowski factor cannot be applied. Figure 2a shows an example of the location of the windings in the transformer window for Yyd 3/6-phase transformers for 12 pulses rectifiers and Figure 2b for Yzyd 3/12-phase transformers for 24 pulses rectifiers. The position of the windings loses the characteristic symmetry for classic three-phase power transformers, see Figure 1. At least the 2D field is necessary to find leakage inductances.
An analytical solution of the magnetic vector potential in transformer air windows with arbitrarily located balanced pairs of windings is presented in [21]. The zero Neumann boundary conditions can be stated for such cases on the ferromagnetic boundaries, which means that the magnetizing current is omitted. We can then calculate the short circuit inductance for any pair of windings. However, this approach cannot be used to determine the leakage inductances of individual windings and the mutual-inductances caused by the stray field. To find these inductances, 2D or 3D field calculations are typically needed, including the entire magnetic circuit with the ferromagnetic core and the air zone. [11][12][13][14][15]22]. This paper presents a 2D numerical approach for determining all leakage inductances σ , L n m , both self, as well as mutual, for arbitrarily located windings considering the magnetic field in the air window only, when surrounded by ferromagnetic walls with high permeability. Difficulties in Rectifier transformers for multi-phase voltages and currents [7] can have windings that are much shorter than the limb's height, and the Rogowski factor cannot be applied. Figure 2a shows an example of the location of the windings in the transformer window for Yyd 3/6-phase transformers for 12 pulses rectifiers and Figure 2b for Y z yd 3/12-phase transformers for 24 pulses rectifiers. The position of the windings loses the characteristic symmetry for classic three-phase power transformers, see Figure 1. At least the 2D field is necessary to find leakage inductances.
An analytical solution of the magnetic vector potential in transformer air windows with arbitrarily located balanced pairs of windings is presented in [21]. The zero Neumann boundary conditions can be stated for such cases on the ferromagnetic boundaries, which means that the magnetizing current is omitted. We can then calculate the short circuit inductance for any pair of windings. However, this approach cannot be used to determine the leakage inductances of individual windings and the mutual-inductances caused by the stray field. To find these inductances, 2D or 3D field calculations are typically needed, including the entire magnetic circuit with the ferromagnetic core and the air zone. [11][12][13][14][15]22,23].
This paper presents a 2D numerical approach for determining all leakage inductances L σ n,m , both self, as well as mutual, for arbitrarily located windings considering the magnetic field in the air window only, when surrounded by ferromagnetic walls with high permeability. Difficulties in formulating the boundary conditions for the case of one winding in the air window has been omitted based on the results of test calculations using the 2D Finite Element Method (FEM). Ferromagnetic plates with sufficiently high magnetic permeability around the air window have been included for that calculations at zero Dirichlet condition on external boundaries. To fulfill the conditions that were followed for those tests, additional ampere-turns have been added at the boundaries. It allowed us to obtain the magnetic vector potential distribution in the air window area the same as form the FEM method, with sufficiently good accuracy. Self-and mutual-leakage inductances can be calculated using an energy approach. This time, however, each pair of windings is represented by two self-inductances and one mutual-inductance instead of one short-circuit inductance. In this way, it is possible to determine all elements of the leakage induction matrix in equation ((1)). When the elementary windings are connected, constituting the final transformer's windings, the matrix of constraints C can be defined and the final inductance matrix takes the form, The paper is organized as follows: In the second section, numerical tests are presented, leading to a formulation of boundary conditions on ferromagnetic walls on the air window for unbalanced ampere-turns in the air window. For this purpose, the 2D FEM has been applied. In Section 3, the finite-difference equations have been formulated for magnetic potential in 2D air zone using the Discrete Partial Differential Operators (DPDO), which are presented in Appendix A. Results have been compared to those from FEM, confirming satisfactory accuracy. Section 4 presents energy-based calculations of self-and mutual-inductances, due to fluxes in the air window. The results for two cases: Windings typical location of a three phase power transformers and for short-circuited part of a winding are presented in Section 5. Finally, conclusions are given in Section 6.

Calculations of Magnetic Field in the Transformer's Air Window with FEM
To determine the leakage inductance matrix L σ , the self-and mutual-inductances for each pair of windings should be calculated. As an example, a pair of arbitrarily located windings is illustrated in Figure 4. formulating the boundary conditions for the case of one winding in the air window has been omitted based on the results of test calculations using the 2D Finite Element Method (FEM). Ferromagnetic plates with sufficiently high magnetic permeability around the air window have been included for that calculations at zero Dirichlet condition on external boundaries. To fulfill the conditions that were followed for those tests, additional ampere-turns have been added at the boundaries. It allowed us to obtain the magnetic vector potential distribution in the air window area the same as form the FEM method, with sufficiently good accuracy. Self-and mutual-leakage inductances can be calculated using an energy approach. This time, however, each pair of windings is represented by two self-inductances and one mutual-inductance instead of one short-circuit inductance. In this way, it is possible to determine all elements of the leakage induction matrix in equation ((1)). When the elementary windings are connected, constituting the final transformer's windings, the matrix of constraints C can be defined and the final inductance matrix takes the form, The paper is organized as follows: In the second section, numerical tests are presented, leading to a formulation of boundary conditions on ferromagnetic walls on the air window for unbalanced ampere-turns in the air window. For this purpose, the 2D FEM has been applied. In Section 3, the finite-difference equations have been formulated for magnetic potential in 2D air zone using the Discrete Partial Differential Operators (DPDO), which are presented in Appendix A. Results have been compared to those from FEM, confirming satisfactory accuracy. Section 4 presents energy-based calculations of self-and mutual-inductances, due to fluxes in the air window. The results for two cases: Windings typical location of a three phase power transformers and for short-circuited part of a winding are presented in Section 5. Finally, conclusions are given in Section 6.

Calculations of Magnetic Field in the Transformer's Air Window with FEM
To determine the leakage inductance matrix σ L , the self-and mutual-inductances for each pair of windings should be calculated. As an example, a pair of arbitrarily located windings is illustrated in Figure 4. Figure 4. 2D problem for a transformer window with two arbitrarily located windings.
Calculation of self-inductances is not a trivial problem because the magnetic field in the air window zone should be generated by one winding only, omitting a magnetic core. It is not possible to formulate the boundary conditions for the equation of the magnetic potential for such a case, as it has been explained in [24]. In 2D, the magnetic potential in the air window is described by the partial differential equation of the form: Calculation of self-inductances is not a trivial problem because the magnetic field in the air window zone should be generated by one winding only, omitting a magnetic core. It is not possible to formulate the boundary conditions for the equation of the magnetic potential for such a case, as it has been explained in [24]. In 2D, the magnetic potential in the air window is described by the partial differential equation of the form: where A(x, y) is the magnetic potential, j 1 (x, y) and j 2 (x, y) are ampere-turns densities of respective windings. The windings with currents i 1 and i 2 have numbers of turns w 1 and w 2 respectively. Taking into account the ferromagnetic, the zero Dirichlet conditions can be determined on the external circumference of the ferromagnetic core. However, by omitting the ferromagnetic core, i.e., assuming µ Fe = ∞, it is not possible to formulate the boundary conditions only for the air window. The zero Neumann conditions are applicable only in the case of balanced ampere-turns To determine the boundary conditions for non-balanced ampere-turns in the transformer window, some numeric test was carried out. The magnetic potential was calculated for the case when only one winding exists in the air window, i.e., w 2 · i 2 = 0, but the ferromagnetic walls are included too. The 2D finite elements software has been used here, assuming zero Dirichlet conditions at the outer boundaries of the ferromagnetic core. The permeability of ferromagnetic µ Fe have been rising up to 10 000 · µ 0 in successive calculations. Figure 5 shows the magnetic potential distribution A(x, y) in the air zone only when one winding generates the field at the highest value of magnetic permeability µ Fe .
x y ∂ ∂ re ( , ) A x y is the magnetic potential, 1 ( , ) j x y and 2 ( , ) j x y are ampere-turns densit ective windings. The windings with currents 1 i and 2 i have numbers of turns 1 w an ectively. Taking into account the ferromagnetic, the zero Dirichlet conditions can be determ e external circumference of the ferromagnetic core. However, by omitting the ferromag , i.e., assuming Fe μ = ∞ , it is not possible to formulate the boundary conditions only for t ow. The zero Neumann conditions are applicable only in the case of balanced ampere [23]. To determine the boundary conditions for non-balanced ampere-turns former window, some numeric test was carried out. The magnetic potential was calculat ase when only one winding exists in the air window, i.e., Energies 2020, 13, 6464 6 of 20 Energies 2020, 13, x FOR PEER REVIEW 6 of 20 The curve ' ' l is the circumference of the air window and ' ' S is the surface limited by this curve. The value of that component equals t In this paper it is proposed to modify the ampere-turns function in the air window instead of solving the Equation (3) That additional surface current is located on the air window boundary, as it is presented in Figure 7. The exciting function is now balanced and satisfies the relation, The curve l is the circumference of the air window and S is the surface limited by this curve. The value of that component equals H t = w 1 · i 1 /(2 · (h + d)) in considered case.
In this paper it is proposed to modify the ampere-turns function in the air window instead of solving the Equation (3) with zero Neumann boundary conditions when omitting j 2 (x, y). This modification results in ampere-turns that are now balanced. To eliminate the ferromagnetic core from the computation a general formula for tangential components on both sides of the boundary H air t − H Fe t = 0 has been changed to H air t − 0 = K z , where K z is a surface current in z direction, which equals to the constant component H t found using FEM That additional surface current is located on the air window boundary, as it is presented in Figure 7.
That additional surface current is located on the air window boundary, as it is presented in Figure 7.   The exciting function is now balanced and satisfies the relation, 0 <x< 2π 0 <y< 2π

Calculations of Magnetic Field in the Transformer Air Window using PDPO
The equation with a modified excitation has the form, where and can be solved for zero Neumann conditions at all boundaries. The modified excitation is compared to the windings ampere-turn function in Figure 8.

Calculations of Magnetic Field in the Transformer Air Window using PDPO
The equation with a modified excitation has the form, where and can be solved for zero Neumann conditions at all boundaries. The modified excitation is compared to the windings ampere-turn function in Figure 8. The solution of that equation can be predicted in the form of the double Fourier series of cosine functions [24], The function mod ( , ) j x y is rather inconvenient for the Fourier analysis, so to solve Equation (6) a special numerical approach has been applied. To find the magnetic potential the DPDO [25,26] have been developed for the series (8). Full consideration leading to the DPDO are presented in Appendix A. The DPDO have been successfully applied for solving 1D boundary problems [24][25][26]. The Finite-Difference Equations (FDE) can be rather easily written using DPDO operators, and it takes the form, The vectors A and J are composed like vectors z and g in Appendix A. The matrix (2) D is singular and has only one eigenvalue equal to zero, so its rank is ((R+1) (S+1) 1) ⋅ − . As it has been mention before the series should satisfy the condition 0,0 A 0 = , i.e., The solution of that equation can be predicted in the form of the double Fourier series of cosine functions [24], The function j mod (x, y) is rather inconvenient for the Fourier analysis, so to solve Equation (6) a special numerical approach has been applied. To find the magnetic potential the DPDO [25,26] have been developed for the series (8). Full consideration leading to the DPDO are presented in Appendix A. The DPDO have been successfully applied for solving 1D boundary problems [24][25][26].
The Finite-Difference Equations (FDE) can be rather easily written using DPDO operators, and it takes the form, The vectors A and J are composed like vectors z and g in Appendix A. The matrix D (2) is singular and has only one eigenvalue equal to zero, so its rank is ((R + 1) · (S + 1) − 1). As it has been mention before the series should satisfy the condition A 0,0 = 0, i.e., It is, in fact, the constraint equation for the elements of the vector A, which can be expressed in the form: where C is the constraint matrix. Finally, the FDEs can be written in the form, The matrix C T · D (2) · C = D (2) C is non-singular, and the new FDEs (12) have a unique solution. The algorithm for solving the FDE (12) has been prepared in MATLAB, and the resulting magnetic potential distribution in the air window area is presented in Figure 9, denoted as DPDO. These results are compared to the results of FEM with a reduced mean value to zero. It is, in fact, the constraint equation for the elements of the vector A , which can be expressed in the form: where C is the constraint matrix. Finally, the FDEs can be written in the form, ( ) The matrix (2) (2) C T ⋅ ⋅ = C D C D is non-singular, and the new FDEs (12) have a unique solution.
The algorithm for solving the FDE (12) has been prepared in MATLAB, and the resulting magnetic potential distribution in the air window area is presented in Figure 9, denoted as DPDO. These results are compared to the results of FEM with a reduced mean value to zero. The magnetic field distribution, has been determined based on ( , ) A x y . The first order DPDO has been used for that.The results, both of FEM and DPDO, are shown in Figure 10.  The magnetic field distribution, H(x, y) = (H x (x, y)) 2 + (H y (x, y)) 2 (13) has been determined based on A(x, y). The first order DPDO has been used for that.The results, both of FEM and DPDO, are shown in Figure 10.
Energies 2020, 13, x FOR PEER REVIEW 8 of 20 It is, in fact, the constraint equation for the elements of the vector A , which can be expressed in the form: where C is the constraint matrix. Finally, the FDEs can be written in the form, ( ) The matrix (2) (2) C T ⋅ ⋅ = C D C D is non-singular, and the new FDEs (12) have a unique solution.
The algorithm for solving the FDE (12) has been prepared in MATLAB, and the resulting magnetic potential distribution in the air window area is presented in Figure 9, denoted as DPDO. These results are compared to the results of FEM with a reduced mean value to zero. The magnetic field distribution, has been determined based on ( , ) A x y . The first order DPDO has been used for that.The results, both of FEM and DPDO, are shown in Figure 10.  The differences between solutions of magnetic field distributions are presented in Figure 11 and seem to be sufficiently small.
Energies 2020, 13, x FOR PEER REVIEW 9 of 20 The differences between solutions of magnetic field distributions are presented in Figure 11 and seem to be sufficiently small. Figure 11. Difference between magnetic field ( , ) H x y resulting from DPDO and FEM. and the solution can be predicted as the Fourier series (7). Application of the DPDO leads to the required solution.

Co-Energy Based Calculations of Self-and Mutual-Leakage Inductances
The inductances  and the solution can be predicted as the Fourier series (7). Application of the DPDO leads to the required solution.

Co-Energy Based Calculations of Self-and Mutual-Leakage Inductances
The inductances L σ n,n , L σ m,m , and L σ n,m in the Formula (1) can be calculated based on field distribution in the air window, obtained with application of the second order DPDOs presented in Appendix A. Two formulae of magnetic co-energy stored in the air can be used here. The first one determines co-energy from the distribution of the magnetic field: (14) and the second is based on lumped inductances of coils, Inductances can be calculated in three steps: -Calculate inductance L σ n,n from the co-energy E co,n at i n = I n and i m = 0, -Calculate inductance L σ m,m from the co-energy E co,m at i n = 0 and i m = I m , -Calculate inductance L σ m,m from the co-energy E co,nm at i n = I n and i m = I m , Energies 2020, 13, 6464 10 of 20 In 2D the formula of the co-energy takes the form, where l eq is an equivalent length of the air window in z-direction, usually it is the circumference taken in the middle of the air window.

Test 1
As the first test, the leakage inductance matrix L σ of a typical power transformer has been determined. The transformer is shown schematically in Figure 1. The inductance matrix can be predicted in the form Design data, more or less related to a power transformer of about 60MVA, 110/60kV, Yd11 have been used for calculations. Its air window and the winding locations are shown in Figure 12  The magnetic vector potential in the air window has been determined by solving equation (12) using the second order DPDO, presented in Appendix A. The modified ampere-turn functions have been determined for individual cases. An exemplary winding '1 is presented in Figure 13. The magnetic field distribution in the air window has been obtained using the first order DPDO. The discretization mesh, with ∆x = 10 mm and ∆y = 50 mm has been chosen, i.e., the number of points in x' direction is 49 and in y direction is 36, which leads to 1764 points in the whole area. The leakage inductances have been obtained from equalities (14)- (16) for an arbitrary value of currents.       The components of magnetic field H x (x, y) and H y (x, y) are shown in Figure 15. The value of H x (x) on the boundaries in direction x is constant and equals H x = ±50A/m and H y (y) on the boundaries in direction y is constant also and equals H y = ±50A/m.
The short circuit inductance, when accounting for the couplings between the windings on the same limb only, equals sc L 82.19 mH = . The same matrix calculated from expressions presented in [24], obtained in 1D approach, has the form, Figure 15. Components H x (x, y) and H y (x, y).
The values of leakage inductances have been calculated using relations (14)- (16) for the number of turns of HV winding '1 . The inductance values are presented in the matrix below. The symmetry of windings has been taken into account, which is reflected in matrix structure.
The short circuit inductance, when accounting for the couplings between the windings on the same limb only, equals L sc = 82.19 mH. The same matrix calculated from expressions presented in [24], obtained in 1D approach, has the form, The corresponding short circuit impedance equals L sc = 86.72 mH. The differences between those matrices seem to be rather high, but the short circuit inductances are almost equal.

Test 2
For the second test, the leakage inductance matrix has been determined when the short circuit of a part of HV winding occurs, as it is presented in Figure 16. The same data as for Test 1 have been assumed, but 20% of HV winding is short circuited. For that test, the other windings have been omitted. For that case, the leakage inductance matrix takes the form, These inductances cannot be found from a 1D approach. To confirm it, the magnetic potential distribution in the short circuit coil case is shown in Figure 17.
assumed, but 20% of HV winding is short circuited. For that test, the other windings have been omitted. For that case, the leakage inductance matrix takes the form, These inductances cannot be found from a 1D approach. To confirm it, the magnetic potential distribution in the short circuit coil case is shown in Figure 17.

Conclusions
This paper presented a method for calculations of self-and mutual-leakage inductances of power transformer windings with arbitrary winding locations in the air windows. It was possible to obtain inductance values based on 2D calculations of the magnetic field only in the air zone, using discrete partial differential operators for 2D periodic functions. The boundary conditions which allowed for the elimination of the transformer ferromagnetic core were determined based on the test calculations of the magnetic field distribution using FEM that included a transformer core of sufficiently high permeability. It followed that the surface current should be located on the boundaries.
To find the solution predicted in the air zone, the discrete partial differential operators of the first-and the second orders have been developed and presented in Appendix A. They allowed us to easily create the finite-difference equations determining the 2D magnetic vector potential distribution in the air window, and consequently, the distribution of the magnetic field strength. The Figure 17. The current density and magnetic potential distribution for short circuited winding '1 .

Conclusions
This paper presented a method for calculations of self-and mutual-leakage inductances of power transformer windings with arbitrary winding locations in the air windows. It was possible to obtain inductance values based on 2D calculations of the magnetic field only in the air zone, using discrete partial differential operators for 2D periodic functions. The boundary conditions which allowed for the elimination of the transformer ferromagnetic core were determined based on the test calculations of the magnetic field distribution using FEM that included a transformer core of sufficiently high permeability. It followed that the surface current should be located on the boundaries.
To find the solution predicted in the air zone, the discrete partial differential operators of the first-and the second orders have been developed and presented in Appendix A. They allowed us to easily create the finite-difference equations determining the 2D magnetic vector potential distribution in the air window, and consequently, the distribution of the magnetic field strength. The results are sufficiently accurate comparing to FEM.
The energy-based method has been successfully used to calculate the leakage inductances, both the self-and mutual-inductances. The matrix of such inductances has been calculated for the design data of a real three-phase power transformer. For a case when a part of one winding on one limb is short circuited, the inductances have also been calculated. Funding: This research, which was carried out under the theme Institute E-2, was funded by the subsidies on science granted by the Polish Ministry of Science and Higher Education.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Discrete Differential Operators for Double Cosine Fourier Series
Discrete Partial Differential Operators (DPDOs) are developed for the two-variable periodic function z(x, y), which is defined in the area 0 < x < 2π, 0 < y < 2π, and can be expanded in the Fourier series with a limited number of terms, z(x, y) = R r=0 S s=0 Z r,s · (cos(r · x/2) · cos(s · y/2)) (A1) Firstly, the unique relations between the function values and the Fourier coefficients can be found. To determine that relation the set of ((R + 1) × (S + 1)) points uniformly distributed in the area has been chosen, as shown in Figure A1.
where ∆x = π/(R + 1),∆y = π/(S + 1). Proper ordering of the point set {x n , y m } allows writing relations between the set of function values {z n,m } and the set of Fourier coefficients {Z r,s } in the form where z is a vector of the function values {z n,m }, ordered, as follows, T wherez n = z n,1 z n,2 · · · z n,S+1  n m , ordered, as follows, The vector Z contains the Fourier coefficients of the series (A1), The T hyper-matrix has the form, The B matrix is square and has dimension ( ) It takes the form, Figure A1. Discretization of the function's area The vector Z contains the Fourier coefficients of the series (A1), The T hyper-matrix has the form, The B matrix is square and has dimension ((S + 1) × (S + 1)). It takes the form, cos(0 · y 1 ) cos(1 · y 1 ) · · · cos(S · y 1 ) cos(0 · y 2 ) cos(1 · y 2 ) · · · cos(S · y 2 ) . . . . . . . . . . . . cos(0 · y S+1 ) cos(1 · y S+1 ) · · · cos(S · y S+1 ) The T matrix has dimension ((R + 1) × (S + 1)). Its inverse matrix can be calculated from the formula, The (T T · T) matrix is a diagonal of the form, The inverse matrix T −1 can be determined from the relation, Obtaining, Energies 2020, 13, 6464 16 of 20 Therefore, the inverse relation to (A3) can be written as, If the first partial derivatives exist, they have the Fourier series forms, (Z x,r,s ) · (sin(r · x) · cos(s · y)) (A7) (Z y,r,s ) · (cos(r · x) · sin(s · y)) (A8) It can be used to find a relation between the Fourier coefficients of the series (A1) and (A7), (A8). They can be written in the forms, Z x and Z y are vectors of the Fourier coefficients of the respective series of partial derivatives (A7), (A8), ordered analogously, where Z x,r = Z x,r,0 Z x,r,1 · · · Z x,r,S and Z y = Z y,0 Z y,1 · · · Z y,R T where Z y,r = Z y,r,0 Z y,r,1 · · · Z y,r,S The matrices R x and R y stand for the differential operators of a 2D periodic function in the frequency domain with respect to 'x' and 'y'. They are diagonal. The matrix R x is constituted by (R + 1) diagonal R x,r matrices, Each of them is diagonal too and has dimension ((S + 1) × (S + 1)) and takes the form, The matrix R y is also constituted by (R + 1) matrices R y,r , R y = diag R y,0 R y,1 · · · R y,R Each matrix R y,r is diagonal too and have (S + 1) elements, R y,r = −diag 0 1 · · · S Analogous to (A3), using the same notations, the relations, Z x = T x · Z x and z y = T y · Z y (A10) can be created. The matrices T x and T y have the forms, sin(0 · x 1 ) · B sin(1 · x 1 ) · B · · · sin(R · x 1 ) · B sin(0 · x 2 ) · B sin(1 · x 2 ) · B · · · sin(R · x 2 ) · B . . . . . . . . . . . .