Magneto-Thermo-Structural Analysis of Power Transformers under Inrush and Short Circuit Conditions

The main equipment responsible for connection and transmission of electric power from generating centers to consumers are power transformers. This type of equipment is subject to various types of faults that can affect its components, in some cases also compromising its operation and, consequently, the electric power supply. Thus, in this paper, electromagnetic, thermal, and structural analysis of power transformers was carried out with the objective of providing the operator with information on the ideal moment for performing predictive maintenance, avoiding unplanned shutdowns. For this, computational simulations were performed using the finite element method (FEM) and, from that, the different transformer operation ways, nominal currents, inrush current, and short-circuit current were analyzed. In this perspective, analyses of the effects that thermal expansion, axial forces, and radial forces exerted were carried out, contributing to possible defects in this type of equipment. As a study object, simulations were carried out on a 50 MVA single-phase transformer. It is important to emphasize that the simulations were validated with real data of measurements and with results presented in the current literature.


Introduction
Considering the great importance of power transformers for energy transmission, several studies have been developed with a focus on the analysis and solution of problems that affect this equipment and put the power system reliability at risk. There are several factors that can generate failures in power transformers during their useful life, but this article will analyze the temperature rise and the deformation of the windings. The increase in the internal temperature in the windings is due to the losses that occur during the operation, most commonly caused by the Joule effect. The deformation of the windings is mostly caused by the action of the magnetic forces generated inside the transformer.
A major aggravating factor in transformer failures is the occurrence of transient phenomena. Among the many transient events that cause failures, those caused by inrush and short-circuit currents are responsible for the majority of operation interruptions in transformers. The magnitude of inrush current may be as high as 10 times or more of transformer-rated current, which can cause malfunction of protection system [1]. Shortcircuit currents induce excessive forces in the transformer windings that result in deformities in the winding, affecting its mechanical and electrical characteristics [2]. They also impose several disturbances such as voltage drop, false disconnections of protection devices, and mechanical stresses in the transformer windings [3].
The correct distinction among the many types of faults is also essential for the correct operation of transformer protection systems. In Yazdani et al. [4], the authors developed The application of FEM makes it possible to analyze complex phenomena in an efficient way, such as in Gołebiowski et al. [15], where the losses due to eddy currents in laminated magnetic nuclei were analyzed. However, despite the fact that FEM presented fairly coherent results, it was necessary to couple different physical phenomena in order to more accurately represent what happened in the equipment, as well as in couplings such as electromagnetic-structural, electromagnetic-thermal, and thermal-fluid, among others.
Another form of coupling is presented in Zhang et al. [16], where a 500 kV 2D finite element elec-tromagnetic-mechanical coupling model was developed to analyze the transient failure process, considering the resistance-temperature relationship and plastic displacement in thermal systems and mechanical winding characteristics.
Several coupling methodologies are already being studied, such as the one presented in Wang et al. [17], a magneto-structural coupling for the short-circuit condition using finite elements, considering the nonlinear characteristics of the spacers in order to assess the distribution of magnetic forces over the windings and deformations generated using a 110 kV transformer as the object of study.
Another coupling methodology was used in Li et al. and Gong et al. [18,19] in order to evaluate the temperature distribution of the transformer winding, since local heating can accelerate the aging of insulating materials. In these works, a magneto-fluidthermal coupling model was developed to calculate the temperature rise and predict the potential critical points of a simplified 2500 kVA dry-type power transformer and a 220 kV transformer.
In Zhang et al. [20], a method was developed to analyze the buckling resistance of transformer windings on the basis of the electromagnetic thermal structure coupling method. The influence of plastic deformation and the thermal effect on the final load of instability is considered in a three-phase dry type transformer, where thermal analysis was developed in 2D and the others in 3D.
Coupling methodologies also involve the interaction between different numerical methods in order to consolidate the results even more precisely, such as the coupling shown in Lima et al. [21], where the trapezoidal method coupled to the FEM was used to analyze the phenomenon of solidary energization in a three-phase bank formed by three single-phase transformers.
Most of the works involve only two physical domains, such as electromagneticthermal, electromagnetic-structural, thermal-fluid, and thermal-structural; the few works that currently combine three physical domains in computer simulations use in at least one of the stages the simulations in 2D. Unlike previous research, in this work, a coupling was made between three different physical domains, magnetic-thermal-structural, where in all stages 3D modeling was used for greater precision and detail of the results.
In this way, using a 50 MVA transformer as a study tool, we developed a 3D magneticthermal-structural analysis via FEM, considering the nominal current, the inrush current, and the short-circuit current to evaluate the heating and the mechanical deformation of the windings. All simulation results were compared with measurement data for validation of the methodology.

Materials and Methods
Before presenting the results obtained in this work and the scientific contributions, we present a synthesis of physical and numerical parameters used to develop this research, depicting some of the main equations used.

Finite Element Method
FEM is a numerical method used to solve complex problems by calculating the nodal energy potential. To apply such methodology, one must divide the complex problem into simpler subproblems ( Figure 1); these simpler subproblems take the form of a geometrical element defined according to the problem to be studied-in the case of this research, the tetrahedral element. One of the main elements used for three-dimensional problems is the tetrahedral element; this element has four nodes, and the point calculation of the nodal potentials can be defined by (1).
From the Cartesian coordinates of the four nodes of each element, one can also calculate the potential according to (2), where D is the determinant of the matrix formed by the coordinates of the element [22].
For the derivation of the shape functions, using the tetrahedral elements, one must use volumetric coordinates, since to model 3D problems, the volume is a quantity to be taken into account; in this way, the Cartesian coordinates can be interpolated through the volume coordinates. To integrate any generic function f(x, y, z) over the entire tetrahedral element, be it in the electromagnetic, thermal, or structural domain, this integration is performed in the reference element according to (3) [23].
where L 1 , L 2 , L 3 are volume coordinates of the element, and J r is the Jacobian matrix needed for the transformation from the volumetric coordinate system to the Cartesian coordinate system.

Magnetic Modeling
The analysis of electromagnetic devices, as in the case of transformers, requires knowledge of electromagnetic phenomena in interior and in the region around the device [24]. To do this, Maxwell's equations are used, which describe the relations between electromagnetic quantities, allowing for temporal and spatial analysis of these types of equipment.
To solve the problem addressed in this work, we applied the magnetic scalar potential to solve (4). B is magnetic flux density and µ is the magnetic permeability of the medium, into formula (5), we obtain div µ grad Ψ = 0 (6) Applying the Galerkin method, we obtain the following relationship [26]: The first integral on the right side of (8) corresponds to the boundary conditions of the problem. Calculating the terms for the discretized domain, This way, the formulation for the problem of the magnetic scalar potential is defined.

Calculation of Losses in Transformer Windings
The total losses P Tot , which results in higher temperatures in transformers, are generated for two kinds of losses that occur during their operation: losses in iron P i and in copper P c . The total loss that occurs in a transformer can be defined from the sum of these distinct losses.
However, for the development of this work, only losses of copper P c , which occurs due to the joule effect, will be considered.
The heat flow static problems address by FEM are heat conduction problems. These problems are represented by a temperature gradient, G, and heat flux density, F [27]. The heat flux density must obey Gauss' Law, which says that the heat flux out of any closed volume is equal to the heat generation within the volume; this law is represented in differential form as div · F = q where q represents volume heat generation. Temperature gradient and heat flux density are also related to one another via the constitutive relationship: where (k) is the thermal conductivity. FEM allows for the variation of conductivity as an arbitrary function of temperature. Usually, the goal is to find the temperature, T, rather than the heat flux density or temperature gradient. Temperature is related to the temperature gradient by Substituting (13) into Gauss' Law and applying the constitutive relationship yields the second-order partial differential equation: FEM solves (14) for temperature T over a user-defined domain with user-defined heat sources and boundary conditions. As in this work the losses by joule effect is modeled as a function of the current density, J, the heat generation can be defined as follows [25]: where σ is the electrical conductivity of the material, ω is the angular frequency, and the current density can be defined taking into account the magnetic vector potential A [28].
The eddy current loss is written in terms of the current density as If A corresponds to the peak value of the vector potential, then the loss in a first order triangle will be where J is the current density, J * is the complex conjugate of the current density, σ is the electrical conductivity of the material, and ∆ is the determinant of the matrix formed by the coordinates of the element. Since A varies linearly over the element, then by (16) so does J; therefore [29], Substituting into (18) Solving the integral, we have This way, it is possible to determine the losses as a function of current density in the coil of transformer.

Calculation of Electromagnetic Forces in Windings
In the structural model, the objective of this work is to analyze the forces acting on transformer windings, which at high values can critically deform the internal structure of windings. According to the electrodynamic theory, force density → f in a given volume of a transformer winding can be defined by the following equation [30,31].
where the magnetic force density   The previous equation can be decomposed into the following relations: The axial force density Radial forces are produced by the magnetic flux density of axial dispersion, presenting different behaviors for the external and internal windings. The value of the axial dispersion magnetic flux density increases from zero in the outermost region of the high voltage winding to a maximum value in the inner diameter (in the interval between the two windings). Therefore, the magnetic flux density of dispersion at the midpoint between the windings can be determined by (24) [32] where B a is the magnetic flux density of axial dispersion, N is the number of turns in the winding, I r is the nominal current of the winding, h w is the height of the winding, and µ 0 is the vacuum permeability. Therefore, the radial force on a medium diameter winding can be determined by (25) as follows [32]: where F r is the total radial force in the winding and D m is the average diameter in the winding.
As for the axial forces, the analytical calculations in a transformer with concentric windings are not simple or precise regarding the calculation of the magnetic flux density of radial dispersion. However, there are methods that provide approximate results, but the axial forces must be analyzed in two distinct ways, designated as ideal and non-ideal conditions [33].
In the ideal condition, the transformers have a uniform distribution of magnetomotive forces in concentric windings of equal length. For this condition, the sum of axial compression close to the midpoint for the external and internal windings can be obtained directly, according to (26) [34,35] where N I r is the magnetomotive force of the windings; D m is the average diameter of the transformer windings, considering the two windings; h w is the height of the windings; d 0 is the space between the windings; and d 1 and d 2 are the radial thicknesses of the external and internal windings, respectively. For the non-ideal condition, there is a significant increase in axial force in these circumstances. Axial forces present greater complexity for the solution through analytical methods. This occurs mainly because of the difficulty in taking into account the curvature of the windings and the presence of the ferromagnetic core [34]. Thus, the density of the mean radial flow in the mean diameter of the transformer is given by (27).
where h e f f is effective length of the radial flow, a is the conductor length, and I max is the nominal winding current.
The axial force on the other winding of the (N I max ) maximum ampere-turns transformer can be determined by means of (28) [36] The action of these forces on the windings in nominal operating conditions already exert a mechanical stress on the windings, and during the occurrence of a transient phenomenon, this stress worsens, further generating deformations in the windings.
The action of radial and axial forces causes different deformations in the windings. The axial force deforms the windings by compressing the upper and lower ends (Figure 3a). The radial force acts more intensely at half the height of the windings, expanding the outer winding and compressing the inner winding (Figure 3b). The mathematical relationship to determine the deformation (in percentage of elongation) of a body in the elastic regime is obtained according to (29) [37].
where ε is the deformation suffered; ∆l is the elongation suffered; and l f − l 0 are the final and initial lengths for the elongation, respectively.

Main Transients in Transformers
Transient phenomena may occur at various times during operation or during energizing of the transformers. These events are configured to raise the peak values of nominal currents in the machine to high levels. The main transients include short-circuit currents and inrush currents. Each phenomenon has specific effects, consequences, and durations as follows.

Inrush
The energizing current has a high initial peak value, as illustrated by Figure 4, and may exceed the peak current of rated current by as much as 20 times. Although this current typically appears during the energization of the transformer, other transients that occur in the circuit may cause this transient current to appear. Therefore, in order to mathematically determine the behavior of the inrush current as a function of time, one can use (30) [38,39]: where φ mp is the maximum magnetic flux; φ r is the residual flow; θ is the angle at the switching instant; N 1 is the number of turns of the winding; and R 1 and L 1 are the resistances and inductances of the primary winding, respectively. However, there are other ways to calculate the peaks of the static inrush current, where the maximum peak is calculated and then coefficients are used to determine the decay of the waveform, in order to determine the current value of the first peak, as in (31): where B mp is the projected steady-state flux density value in the core, B r is the residual flux density, A c is the cross-sectional area of the core, A w is the mean area of a winding loop, and µ 0 is the magnetic permeability of vacuum. In this paper, only the first peak of the inrush condition was analyzed, which is the most critical case for transformers.

Short-Circuit
Short-circuit currents, besides being one of the most frequent causes of faults in transformers, are also among the faults that present greater severity in terms of impact on the structures of transformers support. Under the action of a short-circuit, the dispersion flux density increases significantly and, therefore, the forces acting on the windings also increase [30]. Figure 5 illustrates the characteristic curve of the short-circuit current. For simplification of empirical short-circuit calculations, the impedances of static components such as transmission lines, cables, reactors, and transformers are assumed to be time-invariant. Ignoring these effects and assuming that the transformer impedance (Z) is time-invariant during a short-circuit, the transient and steady-state currents are given by the differential equation of the R − L circuit with an applied sinusoidal voltage [40]: where L is the inductance, R is the resistance, E m is the peak source voltage, ω = 2π f is the angular frequency, f is the frequency of the AC source, and θ is the angle on the voltage wave at which the fault occurs. The solution of this differential equation is given by where I m is the maximum steady-state current, given by E m /Z, and the angle φ = tan −1 (ωL)/R. The transient current, given by the second term of (33), can be called a DC component and it decays at an exponential rate. Equation (33) can be simply written as Although the short-circuit current presents this transient behavior, in this work, only the maximum current point of this I cc phenomenon in a transformer will be analyzed; it can be calculated by (35) [26].
where k is the asymmetry factor of short-circuit current, S n is the nominal power of the transformer in MVA, V is the phase-to-phase nominal voltage of the transformer, and Z is the impedance of the transformer.

Results
The finite element method (FEM) has multiple applications; however, for the present paper, three different types of analysis coupling three different physical domains (magnetic, thermal, and structural) were carried out with three different modes of operation (nominal, inrush, and short-circuit) of a 50 MVA transformer. The simulations were developed using the Ansys Workbench commercial software, a computational tool widely used in academia to develop numerical FEM simulations. The analyses developed in this work are described in the following topics.

Transformer Data
The object of the study was a single-phase 50 MVA transformer connected to a threephase bank that is in operation in Northern Brazil, as illustrated in Figure 6. Its factory characteristics are listed in Table 1.  For the construction of the modeling of the transformer used in this work, we took into account the B-H curve of the ferromagnetic material, where the 1008 cylindrical steel was used, with a packing factor of 0.72. Copper was used to model the windings, in addition to following all the geometric parameters in Table 1.
From the transformer board data, we developed a computational model according to construction parameters. This model is shown in Figure 7.

Simulation in Nominal Operation
In the first moment, the modeling of the transformer under analysis was developed, and then the nominal operating current was inserted as a source of excitation for the electromagnetic analysis. After the solution of the electromagnetic equations, the losses by Joule effect were collected as results. These losses were inserted in a thermal simulation through a coupling, where it was possible to determine the temperature distribution along the transformer windings, as shown in Figure 8. Then, the temperature values obtained were compared with temperature values measured directly inside the transformer at the same time that the current values were collected to validate the developed thermo-magnetic coupling methodology. As shown previously, the losses due to the Joule effect occur as a function of the current density, and thus this type of coupling is necessary to give more precision to the results since the losses are not evenly distributed along the transformer windings.
Considering that the main sources of heat in the transformers are the losses and the Joule effect, we initially performed magnetic simulations by inserting as input data the currents of 131.8 A in the external winding and 421.1 A in the internal winding; these current values were collected directly from the machine in operation. Figure 9 shows the ohmic losses in the transformer windings from these results, and, using coupled simulations, the calculated losses were entered as a boundary condition for the thermal simulation. After conducting the electromagnetic simulation for the nominal operating condition, we found that the losses in the internal winding (Figure 9b) were higher than the losses in the external winding (Figure 9a), thus describing what happens physically in the transformer.
In the sequence, the thermal simulation of the transformer was performed to estimate the temperature of the windings. Ohmic losses calculated in the electromagnetic simulation were inserted in the boundary condition as a heat source; copper emissivity of 0.018, ambient temperature of northern Brazil in the amount of 35 • C, and thermal conversion coefficient were also considered in the amount of 9 W/m 2 • C.
In this way, it was possible to determine the temperatures of the transformer windings for the nominal operating condition, as shown in Figure 10. Analyzing the results of the thermal simulation, we were able to observe that the temperatures in the internal winding ( Figure 10b) were higher than the temperatures in the external winding ( Figure 10a); this was due to the fact that the losses in the internal winding were higher than the losses in the external winding.
In order to validate the computational methodology used in this work, we compared the data obtained with data collected directly from the equipment in nominal operation ( Table 2).  Table 2 shows the summary of the comparison between the thermal data. On the basis of this comparison, one can conclude that the methodology used for the simulations was valid, even having presented a significant margin of error between the values obtained in the simulation and the values measured by temperature sensors in the transformer.
When analyzing the data in Table 2, we were able to observe that the temperature values of the windings measured by sensors were higher than the temperature values found in the simulation. This was due to the fact that in the simulations developed in this work, oil temperatures and losses in the transformer core were not considered, whereas the transformer sensors captured the temperature of the transformer oil and used an analytical methodology to estimate the temperature of the windings.

Simulation with Inrush Current
In the second case, the same modeling of the transformer was used; however, in this case, the maximum peak value of the inrush current was applied as the energizing current for the electromagnetic simulation. After processing, the axial and radial forces generated in the external winding were collected as a result, as shown in Figure 11. Then, the calculated forces were inserted as a boundary condition in the structural simulation through the magneto-structural coupling, where it was possible to determine the radial deformation in the external winding. The deformation results found were compared with results found in the literature in order to validate the magneto-structural coupling methodology developed in this work. Once the thermal part methodology was validated, simulations were performed, considering more severe operating conditions in order to analyze the deformation caused by the action of magnetic forces inside the machine. Figure 12 shows the results obtained in magnetic simulation for inrush current operation, in which case a current of 1733.1 A was used in the external winding. From the axial and radial magnetic flux density, we were able to calculate the radial and axial forces, respectively. The calculated forces were entered as a boundary condition for the structural simulation in order to estimate winding deformation, as shown in Figure 13. For the inrush condition, only the external winding was energized, and thus there was only the generation of electromagnetic forces in one of the windings.
Observing the results of the structural simulation, we obtained the value of 6.4 × 10 −3 mm as the maximum radial strain. The values found in this stage of the simulations were higher than the values found in the literature, as for example in Fonseca et al. [13], where a similar methodology was used with the same equipment under the same operating conditions. However, in Fonseca et al. [13], the electromagnetic stage was modeled in 2D, while in this work the electromagnetic part was modeled in 3D.
The difference between the results can be justified in Feng et al. [41], where the author compared simulation computational results in 2D and 3D and showed that, even with acceptable values, 3D simulations always present values higher than 2D results. This is because in 2D models, a simplification is used in geometries that present an axis of symmetry in addition to the use of a depth coefficient to approximate the model, whereas in 3D cases, the problem is fully modeled in the three dimensions.
Therefore, the methodology for the simulations can be considered valid.

Simulation with Short-Circuit Current
One of the main contributions developed in this article that is characterized by performing a 3D multiphysics analysis of a power transformer for the short-circuit condition coupling electromagnetic-thermal-structural phenomena. After validating the previous coupling methodologies, we built a new methodology based on the previous oneshowever, in this case, using the maximum peak current of the short-circuit current as the excitation source.
Using the same modeling of the transformer shown above and inserting the shortcircuit current as a source of excitation for the electromagnetic simulation, after processing, we extracted the Joule losses and axial and radial forces as results. In the sequence, the losses calculated in the electromagnetic analysis were inserted in the thermal simulation as a boundary condition, where it was possible to determine the temperature distribution along the windings of the transformer, and it was also possible to determine the thermal expansion suffered by the windings due to the high temperatures. Finally, a structural analysis was coupled to the set where the axial and radial forces calculated in the elec-tromagnetic analysis and the thermal expansion calculated in the thermal analysis were inserted as a boundary condition. In this way, it was possible to reach a more precise structural deformation of the windings to the transformer since these consider electromagnetic and thermal elements that directly influence the structure of the transformer. The developed scheme for coupling the electromagnetic, thermal, and structural physical domains is illustrated in Figure 14. As the last condition for the simulations of this work, short-circuit currents were adopted, a requirement that is considered to be one of the most severe among transients. Under these conditions, currents with amplitudes of 8979 A were adopted in the external winding and 34,586 A in the internal winding, where the following results were obtained.
Like the previous step, in the electromagnetic simulation, the current density was calculated in windings during the maximum peak of a short-circuit; due to Joule losses, these high-current values flowing through a conductor cause extreme heating inside the equipment, as shown in Figure 15. Analyzing the temperature shown in Figure 15, we were able to see how critical a short-circuit is to a transformer, reaching temperatures that can break insulation or even melt winding discs. At this point, it is important to note that the current protection system should quickly isolate the fault, preventing the equipment's temperature from reaching these extremely high levels. In this paper, this protection was not considered exactly for demonstrating the temperature level that can be reached if the protection system fails, which would certainly lead the transformer to irreparable damage.
In addition to current density, the magnetic simulations provide the determination of magnetic forces generated inside the transformer during short-circuit. From the magnetic forces and using coupled simulations, we determined the deformation suffered due to the interaction of magnetic forces and thermal expansion due to the heating. The joint action of all these factors resulted in the deformation illustrated in Figure 16. According to Figure 16, it is possible to observe that the deformation was more critical in the external winding; considering that the two previous conditions were validated, the results obtained in this step can be considered valid. Although the maximum deformation achieved reached 2.15 mm, this would already be enough to break the insulation and drastically compromise the operation of the transformer.

Conclusions
In this work, several analyses were developed combining electromagnetic, thermal, and structural phenomena using the FEM to analyze the behavior of a 50 MVA transformer operating with nominal current, inrush current, and short-circuit current.

1.
Initially the transformer was simulated using a thermo-magnetic methodology and considering the nominal operating current. The results obtained were compared with values measured by sensors to validate the methodology.

2.
In the second case, a magneto-structural methodology was used to simulate the same equipment; however, in this case, the maximum peak of the inrush current was considered as a source of excitation. Even though the results obtained in this step were higher than the results found in the literature, this difference is perfectly justified in view of the way in which this methodology was constructed.

3.
In the last step of the methodology of this work, a magneto-thermo-structural methodology was developed where the 50 MVA transformer was simulated considering the maximum peak of the short-circuit current as an operating condition. In this case, it was considered the most severe case that can occur in the transformers. The results found showed internal temperatures that can melt the contours, in addition to deformations that can easily break the insulation of the turns to the transformer. 4.
The methodology used in this work is very promising and the results obtained describe a projection for the physical behavior of the transformer; however, there are still points to be improved in future works, such as the use of this coupling methodology in transient analyses considering the operation curves, as well as the variation of the phenomena over time. Informed Consent Statement: Not applicable.

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

Nomenclature
FEM Finite element method HTS High-temperature superconducting SVP Saturated vapor pressure CTC Continuously transposed conductors