Numerical Investigations for Vibration and Deformation of Power Transformer Windings under Short-Circuit Condition

: The analysis of the dynamic process of winding destabilization under sudden short-circuit conditions is of great importance to accurately assess the short-circuit resistance of power transformers. Based on magneto-solid coupling, an axisymmetric model of the transformer and a 3D multilayer model of the transformer considering the support components are established, respectively, and the short-circuit electromagnetic force (EF) is simulated by using the ﬁnite element method. It is concluded that the middle layer of the winding is subjected to the larger radial EF, while the axial EF has a greater effect on the layers at both ends. Moreover, the impression of the preload force, aging temperatures, and the area share of spacers on the vibration and deformation of windings are studied under short-circuit conditions. The overall distribution of plastic strain and residual stress in the winding is symmetrical, and the maximum values occur in the lower region of the middle of the winding. Finally, considering the material properties of disks and insulating components, the cumulative effect of plastic deformation under multiple successive short-circuit shocks is calculated. Compared with the traditional axisymmetric model of transformer, the three-dimensional multilayer model of the transformer established in this paper is more suitable for the actual winding structure and the obtained results are more accurate.


Introduction
With the improvement in the voltage level of power systems, short-circuit accidents caused by the damage of the transformer have been becoming more and more serious [1,2].The stability and safety of power transformers determines the reliability of the power system.When a short-circuit fault occurs, the short-circuit current is approximately dozens of times the rated current, which leads to a short-circuit electromagnetic force (EF) hundreds of times stronger than that under normal operating conditions [3].Under such strong EF, transformer windings are highly susceptible to accidents such as deformation, inter-turn insulation breakage, and instability.
We refer to the coupling between the EF on the winding and its structural deformation as magneto-mechanical coupling.A considerable number of researchers have developed transformer models using the finite element method to study the magneto-structure coupling effect of windings [4][5][6][7].The researchers equated transformer disks as concentrated masses and insulation spacers as springs to study the axial vibration of the windings based on the mass-spring model [8,9].After a short-circuit shock in the transformer, the winding deformation affects both the current distribution and the leakage flux density distribution, thus changing the EF acting on the winding, which in turn affects the winding deformation.Recently, many research results have also shown that the deformation and structure of windings have a non-negligible effect on EF [10][11][12][13].
The main factors affecting the vibration and deformation of transformer windings include the nonlinearity of the insulation spacer material [14] and the axial preload force.In Energies 2023, 16, 5318 2 of 18 our work, magneto-mechanical coupling models of the power transformer are established to investigate the distribution of leakage flux and short-circuit EF as well as factors of the winding vibration.In addition, the cumulative effect of plastic deformation under multiple successive short-circuit shocks is analyzed.

Electromagnetic Analysis
The completion of electromagnetic analysis is twofold.Firstly, the short-circuit current is extracted based on a circuit model of transformers.Secondly, the short-circuit current is imposed to the windings as the excitation of Maxwell's equation to calculate the detailed distributions of both leak flux and electromagnetic force.
It is worth emphasizing that when short-circuit faults occur, the cores of the transformers are not saturated [15]; thus, we can assume a constant short-circuit impedance and start with the equivalent model shown in Figure 1, which is governed by where u 1 is the primary voltage, U 1 is the root mean square (RMS) of u 1 , ω is the angular frequency, t is the time, α is the initial phase angle, i k is the transient short-circuit current, and R k and L k are the short-circuit resistance and inductance, respectively.
Energies 2023, 16, x FOR PEER REVIEW 2 of 19 The main factors affecting the vibration and deformation of transformer windings include the nonlinearity of the insulation spacer material [14] and the axial preload force.In our work, magneto-mechanical coupling models of the power transformer are established to investigate the distribution of leakage flux and short-circuit EF as well as factors of the winding vibration.In addition, the cumulative effect of plastic deformation under multiple successive short-circuit shocks is analyzed.

Electromagnetic Analysis
The completion of electromagnetic analysis is twofold.Firstly, the short-circuit current is extracted based on a circuit model of transformers.Secondly, the short-circuit current is imposed to the windings as the excitation of Maxwell's equation to calculate the detailed distributions of both leak flux and electromagnetic force.
It is worth emphasizing that when short-circuit faults occur, the cores of the transformers are not saturated [15]; thus, we can assume a constant short-circuit impedance and start with the equivalent model shown in Figure 1, which is governed by where u is the primary voltage,  The transient short-circuit current () k it can then be calculated as where k I is the RMS of steady-state current.Once () k it is obtained, it can be used to excite the Magnetoquasistatic (MQS) field governed by The transient short-circuit current i k (t) can then be calculated as where I k is the RMS of steady-state current.Once i k (t) is obtained, it can be used to excite the Magnetoquasistatic (MQS) field governed by as well as the constitutive law Energies 2023, 16, 5318 3 of 18 In the above equations, A is the magnetic vector potential, E is the electric field intensity, H is the magnetic field intensity, B is the magnetic flux density, J is the current density, J e is the applied current density, v is the velocity of the conductor, σ is the electrical conductivity, and µ is the magnetic permeability.A pregiven J e is imposed on the ports of windings.
After the extraction of the leakage flux via solving Equations ( 3) and ( 4), the resultant volume density of the short-circuit EF is given by (5)

Mechanical Analysis
For accurate analysis of the elastic-plastic deformation and the cumulative effect of transformer winding, it is necessary to determine the internal stress state induced by EF on the windings.The short-circuit EF is a time-varying and uneven physical force exciting the transient structural equation where is the stress tensor ρ is the mass density, u = [u x , u y , u z ] T is the displacement vector.It is worth noting that in Equation ( 7) τ ij = τ ji , where i, j ∈ {x, y, z}.Similarly, the strain tensor is also a symmetric tensor.The two tensors are linked through the constitutive laws where E and v denote the elasticity modulus and Poisson's ratio, respectively.The displacement vector u is related to ε by the geometric law

Finite Element Solution to the Coupled Magneto-Mechanical System
The finite element method (FEM) is adopted for the numerical solution to the coupled magneto-mechanical equations due to its capability of resolving complex geometry.The resultant discretized system in a compact matrix form can be described by where M is the mass matrix, C is the damping matrix, K(U) is the stiffness matrix depending on the displacement vector U, F(A) is the EF force vector calculated according to the magnetic vector potential A, and K m is the coefficient matrix representing the MQS equation.Note that the dependence of K(U) on U makes Equation (12) a nonlinear system.The nonlinearity stems from nonlinear mechanical characteristics of some specific structures, e.g., the insulating spacers.This nonlinearity will be defined in detail in Section 3.
Another point is that although Equation (3) involves the first-order time derivative of E, the discretized MQS equation does not contain any time derivative, because we choose to solve only the magnetic vector potential A, which behaves like a static field under the MQS assumption.
The FEM discretization of the governing equations is realized with the commercial FEM solver COMSOL Multiphysics 6.0, and we choose the built-in backward Euler algorithm allowing adaptive time step sizes to complete the time advances of Equation (12).
As for the implementation details of the FEM discretization, including the shape functions, evaluation of the coefficient matrices, and time-stepping algorithms, the reader can refer to [16].

Description of the Model
In this section, the vibration and deformation of the windings under the short-circuit condition of a three-phase oil-immersed-type 110 kV power transformer are thoroughly investigated as an example.The transformer has three groups of windings, the rated voltages of which are 110 kV, 38.5 kV, and 10.5 kV.In this work, we focus on the interactions of the latter two groups of windings, referred to as the HV winding and LV winding, respectively.The HV windings are under a star arrangement, while the LV windings are delta-connected.Some key parameters of the transformer are given in Table 1.An axisymmetric model of the transformer based on the spring-mass system and a 3D multilayer model of the transformer considering the support components are established, as shown in Figure 2, respectively.The force between the adjacent windings of the threephase transformer is neglected under a sufficient axial preload, thus we only calculate one phase winding [17].The electromagnetic field and mechanics share the same set of meshes and have an identical form of mesh division.There are, of course, some similarities and differences between the two finite element models in Figure 1.
In Figure 2a, a simplified axisymmetric model of the transformer is presented, which is obtained by neglecting some unsymmetric geometry details.The axisymmetric model is favorable regarding computational efficiency and suitable for the simulation of the axial vibration of the winding.In Figure 2b, the 3D multilayer model reflects the actual spacers' distribution and fits better to the real winding structure.Meanwhile, the number of turns Energies 2023, 16, 5318 5 of 18 of primary and secondary windings remains the same with the parameters in Table 1, and the number of disks is simplified to some extent to save calculation time.
An axisymmetric model of the transformer based on the spring-mass system and 3D multilayer model of the transformer considering the support components are estab lished, as shown in Figure 2, respectively.The force between the adjacent windings of th three-phase transformer is neglected under a sufficient axial preload, thus we only calcu late one phase winding [17].The electromagnetic field and mechanics share the same se of meshes and have an identical form of mesh division.There are, of course, some simi larities and differences between the two finite element models in Figure 1.As for material properties, the core is made of soft iron with a nonlinear B-H curve, which is depicted in Figure 3.The windings are made of copper with σ = 5.7 × 10 7 S/m, µ r = 1, E = 110 × 10 9 Pa, and v = 0.34.For disks and insulating spacers, we set σ = 1 S/m, µ r = 1, and v = 0.30.The nonlinear elastic modulus is specified in Section 3.4.1.

Calculation of the Leakage Flux and Short-Circuit EF
Here, the short-circuit current is considered when a three-phase symmetric shortcircuit fault occurs.After a short circuit occurs, the peak short-circuit current reaches its maximum, which is approximately 7 × 10 4 A, at t = 0.01 s.The transient short-circuit current obtained by Equation ( 2) is shown in Figure 4, and it is applied to the winding as an excitation source.
The leakage flux at the time instant t = 0.01 s after the sudden short circuit, i.e., when the instantaneous value of the short-circuit current reaches its maximum, is shown in Figure 5.The short-circuit EF is affected by both the short-circuit current and the leakage field.The distributions of EF in HV and LV windings at 0.01 s are given in Figure 6.It is worth noting that in our 3D simulation, the number of tetrahedral elements is 164443, which is sufficient to yield almost grid-independent solutions for our example.Therefore, the simulation results here and in the following sections are reliable from the numerical perspective.the number of disks is simplified to some extent to save calculation time.
As for material properties, the core is made of soft iron with a nonlinear Bwhich is depicted in Figure 3.The windings are made of copper with 5.7

Calculation of the Leakage Flux and Short-Circuit EF
Here, the short-circuit current is considered when a three-phase symmetr circuit fault occurs.After a short circuit occurs, the peak short-circuit current re maximum, which is approximately   field.The distributions of EF in HV and LV windings at 0.01 s are given in Figure worth noting that in our 3D simulation, the number of tetrahedral elements is 1 which is sufficient to yield almost grid-independent solutions for our example.The the simulation results here and in the following sections are reliable from the num perspective.The overall radial component of winding EF is strong in the middle and mild at both ends, due to the concentrated distribution of axial leakage in the middle.The negative radial component of the short-circuit EF in the primary winding means that the primary winding is compressed inward, while the positive radial component of the short-circuit EF in the secondary winding indicates that the secondary winding is expanded outward.As the radial leakage is concentrated at both ends, the axial EF is at its maximum at the end of the winding.In addition, the closer it is to the middle of the winding, the smaller the axial force is.

Effect of the Area Proportion of Spacers on Axial Vibration
The contact area between the disk and the spacer is an essential indicator for the design of winding short-circuit resistance.Because all insulation spacer layers have the same structure, only one of them is illustrated.Figure 7 shows a schematic cross-section of a particular insulating spacer layer along the axial direction, where the white area is the air and the gray area indicates the insulating pads.By varying the width of each spacer to change the contact area between each layer of the disks and the insulating spacers, the study index is set to the percentage ϕ 1 of the spacer contact area over the circumference of the spacers.Considering the insulating spacer layer as a composite material consisting of a insulating paperboard, the composite elasticity modulus can be expressed as: where  1 is the area proportion of spacers in each layer, E 1 is the elasticity modu the cushion block material,  2 is the area proportion of air in each layer, and E 2 elasticity modulus of air.If the contact areas of the spacers and the wire cakes va ,  2 , and the resultant equivalent elasticity modulus of the winding layers will c accordingly.
In Figure 8, the horizontal coordinates indicate the numbering of the studied The neighboring biscuits are numbered in order from top to bottom.For instance indicates that the studied objects are the two biscuits adjacent to the uppermost end winding, and No.2 represents that the studied objects are the second and third b counted from the upper end.The LV winding has a total of 76 bobbins, so there numberings in total.The HV winding has a total of 86 bobbins, resulting in a tota numbers.Considering the insulating spacer layer as a composite material consisting of air and insulating paperboard, the composite elasticity modulus can be expressed as: where ϕ 1 is the area proportion of spacers in each layer, E 1 is the elasticity modulus of the cushion block material, ϕ 2 is the area proportion of air in each layer, and E 2 is the elasticity modulus of air.If the contact areas of the spacers and the wire cakes vary, ϕ 1 , ϕ 2 , and the resultant equivalent elasticity modulus of the winding layers will change accordingly.
In Figure 8, the horizontal coordinates indicate the numbering of the studied object.The neighboring biscuits are numbered in order from top to bottom.For instance, No.1 indicates that the studied objects are the two biscuits adjacent to the uppermost end of the winding, and No.2 represents that the studied objects are the second and third biscuits counted from the upper end.The LV winding has a total of 76 bobbins, so there are 75 numberings in total.The HV winding has a total of 86 bobbins, resulting in a total of 85 numbers.
indicates that the studied objects are the two biscuits adjacent to the uppermost end winding, and No.2 represents that the studied objects are the second and third b counted from the end.The LV winding has a total of 76 bobbins, so there numberings in total.The HV winding has a total of 86 bobbins, resulting in a tota numbers.The vertical coordinate represents the difference between the displacements two adjacent spacers.More specifically, the longitudinal coordinate indicates the ence in displacement between the lower plane of the upper one and the upper p the lower one in the adjacent disks.The ordinate coordinate is denoted by D .If there is no relative displacement between the reference faces representing the ad two disks, and the distance between them is the width of the spacer.If D 0  , t tance between the reference faces representing adjacent disks is less than the width spacer.That is, the spacer is compressed, the disk is in close contact with the space there is no gap between them.When D 0  , the distance between the reference su representing adjacent disks is greater than the width of the spacer, i.e., a gap is c between the disks and the spacer, which has the risk of winding instability.
From Figure 8a, the longitudinal coordinates are all less than zero, which mea the disks and spacers are in close contact, and the stability of the LV windings is sa tory.From Figure 8b, the gaps appear between the top and bottom disks and spa HV windings.The end disks are prone to instability, while the middle disks hav stability.Moreover, the stability of both LV and HV windings changes most wh The vertical coordinate represents the difference between the displacements of the two adjacent spacers.More specifically, the longitudinal coordinate indicates the difference in displacement between the lower plane of the upper one and the upper plane of the lower one in the adjacent disks.The ordinate coordinate is denoted by D. If D = 0, there is no relative displacement between the reference faces representing the adjacent two disks, and the distance between them is the width of the spacer.If D < 0, the distance between the reference faces representing adjacent disks is less than the width of the spacer.That is, the spacer is compressed, the disk is in close contact with the spacer, and there is no gap between them.When D > 0, the distance between the reference surfaces representing adjacent disks is greater than the width of the spacer, i.e., a gap is created between the disks and the spacer, which has the risk of winding instability.
From Figure 8a, the longitudinal coordinates are all less than zero, which means that the disks and spacers are in close contact, and the stability of the LV windings is satisfactory.From Figure 8b, the gaps appear between the top and bottom disks and spacers of HV windings.The end disks are prone to instability, while the middle disks have good stability.Moreover, the stability of both LV and HV windings changes most when ϕ 1 varies from Energies 2023, 16, 5318 10 of 18 20% to 40%.In summary, the area proportion of spacers ϕ 1 boosts the stability of both LV and HV windings; the vibration stability of LV windings is better than that of HV windings.
In the design of the short-circuit resistance of the winding, an appropriate ϕ 1 is essential to avoid the situation of winding instability.Moreover, the fatigue life of the insulation layers of the windings near the ends under continuing transient shock loads should be properly strengthened, so that the breakage due to the frequent collision between the end disks and the spacers can be reduced.Possible measures may be adjustments to the material or thickness.

Effect of Axial Preload on Axial Vibration
The axial preload is an essential parameter in the study of winding short-circuit processes.The effect of axial preload on the axial stability of disks is mainly reflected in changing the elasticity modulus of spacers.The simulation results, obtained by changing the axial preload applied to both ends of the winding for the same winding structure, are shown in Figure 9.In Figure 9, the vibration stability of LV and HV windings improved with the g of the axial preload from 3 MPa to 6 MPa.LV windings are less prone to instability the axial preload is greater than the value between 3 and 4 MPa, while HV windin In Figure 9, the vibration stability of LV and HV windings improved with the growth of the axial preload from 3 MPa to 6 MPa.LV windings are less prone to instability when the axial preload is greater than the value between 3 and 4 MPa, while HV windings undergo stable vibration when the axial preload is greater than 5 MPa.However, high axial preload may cause the windings to tilt or even collapse.
For the transformer model used in this paper, the axial preload should be deliberately selected to ensure the winding axial vibration stability.As seen from Figure 10, the axial displacement differences of the primary an ondary windings of each disk remain negative at all temperatures.It is indicated t distance between the two disks at each moment in the vibration process is less th As seen from Figure 10, the axial displacement differences of the primary and secondary windings of each disk remain negative at all temperatures.It is indicated that the distance between the two disks at each moment in the vibration process is less than the thickness of the insulation spacers.That is, the insulation spacer is always in a compressed state, and the middle spacers of the windings are compressed significantly more than the spacers at both ends of the winding.There is no gap generated between each disk and insulation spacer, and the axial vibration of transformer windings is a stable dynamic process during the short circuit at the above temperature conditions.Overall, the Young's modulus value of the insulation material decreases gradually with the increase in aging temperature, and the degree of the pad is gradually compressed and reduced, which demonstrates that the axial stability of the winding degrades gradually.The simulation results are verified against the experimental conclusions of the existing research.
Using the 3D multilayer model shown in Figure 2b, the winding residual stress distributions under different aging temperatures are obtained when the short-circuit EF amplitude reaches the first wave peak, as shown in Figure 11.In Figure 11, at the maximum amplitude of the short-circuit EF, the stress value of the disks in the middle of the winding is greater than that of the end disks under the radial short-circuit EF.The maximum surface stresses of the primary and secondary windings at 20 °C, 60 °C, and 130 °C are 1.698 × 10 8 N/m 2 , 1.499 × 10 8 N/m 2 , and 1.116 × 10 8 N/m 2 , respectively, and they are all distributed on the primary winding near the support parts.As the temperature increases, the Young's modulus of the insulating paperboard decreases, the elastic support stiffness of the support components to the winding diminishes, and the maximum stress value on the surface of the winding drops gradually.The key to the determination of the spoke destabilization lies in the critical stress value which is equivalent to the same short-circuit EF.The maximum stress value that the winding can withstand decreases gradually as the aging temperature rises, and the transformer is more prone to radial destabilization accidents.In Figure 11, at the maximum amplitude of the short-circuit EF, the stress value of the disks in the middle of the winding is greater than that of the end disks under the radial short-circuit EF.The maximum surface stresses of the primary and secondary windings at 20 • C, 60 • C, and 130 • C are 1.698 × 10 8 N/m 2 , 1.499 × 10 8 N/m 2 , and 1.116 × 10 8 N/m 2 , respectively, and they are all distributed on the primary winding near the support parts.As the temperature increases, the Young's modulus of the insulating paperboard decreases, the elastic support stiffness of the support components to the winding diminishes, and the maximum stress value on the surface of the winding drops gradually.The key to the determination of the spoke destabilization lies in the critical stress value which is equivalent to the same short-circuit EF.The maximum stress value that the winding can withstand decreases gradually as the aging temperature rises, and the transformer is more prone to radial destabilization accidents.

Mechanical Properties of Disks and Spacers
Once the internal mechanical stress value of the windings exceeds the yield strength of the copper conductor, the windings of the transformer produce a slight plastic deformation in the process of the short-circuit fault.In terms of multiple continuous short-circuit currents, eventually, the gradually accumulated plastic deformation leads to winding deformation, dislocation, and other accidents.
The stress-strain curve of the copper conductor in the transformer winding, as a plastic material, is shown in Figure 12.When a short-circuit fault occurs, the copper conductor enters the elastic deformation stage during the destruction of the winding and begins to produce plastic strain after reaching the yield limit, and finally, the plastic deformation accumulates gradually under multiple short-circuit force impacts.The spacers are made of nonlinear materials, which means that the depen the stress  on the strain  exhibits a nonlinear behavior.For the insulatin the  -expression within a stress of 1000 MPa can be expressed by Equation the variation of the elasticity modulus with strain is shown in Equation ( 15) [18] 3 a b where a is a constant and b is the hardening coefficient.Typically, a is 105 b is 1750 MPa.It is clear from Equation ( 15) that the elasticity modulus E of t The spacers are made of nonlinear materials, which means that the dependence of the stress σ on the strain ε exhibits a nonlinear behavior.For the insulating spacer, the σ − ε expression within a stress of 1000 MPa can be expressed by Equation ( 14), and the variation of the elasticity modulus with strain is shown in Equation ( 15) [18].
where a is a constant and b is the hardening coefficient.Typically, a is 105 MPa and b is 1750 MPa.It is clear from Equation ( 15) that the elasticity modulus E of this material depends on strain.In addition, the relationship between the elasticity modulus and strain can be visualized in Figure 13.Equivalent plastic strain and residual stress are used to represent the strain and stress inside the windings, and a single variable is used to express the degree of plastic deformation.Once the plastic strain is generated, residual stresses are caused by the uneven plastic deformation inside the windings.The residual stress is the internal stress, that remains in the windings after the unloading of the external short-circuit electromagnetism force, and it is in self-balance.When the residual stress and the mechanical stress caused by the external forces superimposed on each other exceed the strength limit of the material, it is likely to cause warping, cracking, or twisting deformations of the structure.If the loading and unloading process of short-circuit EF inside the lined cake generates a different degree of residual stress, it will seriously endanger the stable operation of transformer.

Cumulative Effect of Plastic Strain under Single Short-Circuit Shock
The duration of a short-circuit shock is from t = 0 s to t = 0.1 s after the short circuit occurs.During a single short-circuit shock, each wave of short-circuit EF may cause plastic deformation of the copper winding.The calculated results of plastic strain and residual stress are shown in Figures 14 and 15 for multiple wave-fronts under a single short-circuit shock, respectively, mainly including the t = 0.01 s instant when the short-circuit EF amplitude reaches its maximum and the t = 0.1 s instant after the end of the single shortcircuit shock.Equivalent plastic strain and residual stress are used to represent the strain and stress inside the windings, and a single variable is used to express the degree of plastic deformation.Once the plastic strain is generated, residual stresses are caused by the uneven plastic deformation inside the windings.The residual stress is the internal stress, that remains in the windings after the unloading of the external short-circuit electromagnetism force, and it is in self-balance.When the residual stress and the mechanical stress caused by the external forces superimposed on each other exceed the strength limit of the material, it is likely to cause warping, cracking, or twisting deformations of the structure.If the loading and unloading process of short-circuit EF inside the lined cake generates a different degree of residual stress, it will seriously endanger the stable operation of transformer.

Cumulative Effect of Plastic Strain under Single Short-Circuit Shock
The duration of a short-circuit shock is from t = 0 s to t = 0.1 s after the short circuit occurs.During a single short-circuit shock, each wave of short-circuit EF may cause plastic deformation of the copper winding.The calculated results of plastic strain and residual stress are shown in Figures 14 and 15 for multiple wave-fronts under a single short-circuit shock, respectively, mainly including the t = 0.01 s instant when the shortcircuit EF amplitude reaches its maximum and the t = 0.1 s instant after the end of the single short-circuit shock.As can be seen in Figure 14, the maximum magnitudes of the short-circuit current and short-circuit EF are reached at t = 0.01 s, corresponding to a maximum plastic strain of 0.449%.The maximum value of plastic strain increases to 1.571% at t = 0.1 s after the end of the single short-circuit shock.The gradual increase in the plastic strain maximum with time shows that multiple wave crests under a single short-circuit shock have a certain cumulative effect.Meanwhile, the results at different time instants indicate that the distribution of plastic strain tends to become more uniform as time advances.From Figure 15, the maximum value of residual stress increased from 7.446 × 10 7 N/m 2 to 8.480 × 10 7 N/m 2 during the duration of a single short-circuit shock, and they all appeared on the disks of the primary winding near the end.The distribution of residual stresses is similar to the distribution of magnetic field leakage under short-circuit conditions as well as the distribution of short-circuit EF.
In the radial direction, the residual stress is concentrated between the two windings near the gap.In the axial direction, it is concentrated in the middle of the windings.The primary winding is subjected to a larger radial short-circuit EF than the secondary winding, so the overall residual stress value is higher than that of the residual stress of the secondary winding.In addition, the minimum value of residual stresses exists and is not zero because of the restraining effect of preload and gravity on the winding in the steady state.

Cumulative Effect of Plastic Strain under Multiple Consecutive Short-Circuit Shocks
Due to the multiple reclosures of power systems, recurring short-circuit shocks are a realistic risk to transformers.In this subsection, the stress-strain results obtained from a As can be seen in Figure 14, the maximum magnitudes of the short-circuit and short-circuit EF are reached at t = 0.01 s, corresponding to a maximum plast of 0.449%.The maximum value of plastic strain increases to 1.571% at t = 0.1 s end of the single short-circuit shock.The gradual increase in the plastic strain m with time shows that multiple wave crests under a single short-circuit shock have cumulative effect.Meanwhile, the results at different time instants indicate that th bution of plastic strain tends to become more uniform as time advances.From Figure 15, the maximum value of residual stress increased from 7.4 N/m 2 to 8.480 × 10 7 N/m 2 during the duration of a single short-circuit shock, and appeared on the disks of the primary winding near the end.The distribution of stresses is similar to the distribution of magnetic field leakage under short-circu tions as well as the distribution of short-circuit EF.
In the radial direction, the residual stress is concentrated between the two w near the gap.In the axial direction, it is concentrated in the middle of the windin primary winding is subjected to a larger radial short-circuit EF than the secondar ing, so the overall residual stress value is higher than that of the residual stres secondary winding.In addition, the minimum value of residual stresses exists an zero because of the restraining effect of preload and gravity on the winding in th state.

Cumulative Effect of Plastic Strain under Multiple Consecutive Short-Circu Shocks
Due to the multiple reclosures of power systems, recurring short-circuit sho realistic risk to transformers.In this subsection, the stress-strain results obtained As can be seen in Figure 14, the maximum magnitudes of the short-circuit current and short-circuit EF are reached at t = 0.01 s, corresponding to a maximum plastic strain of 0.449%.The maximum value of plastic strain increases to 1.571% at t = 0.1 s after the end of the single short-circuit shock.The gradual increase in the plastic strain maximum with time shows that multiple wave crests under a single short-circuit shock have a certain cumulative effect.Meanwhile, the results at different time instants indicate that the distribution of plastic strain tends to become more uniform as time advances.
From Figure 15, the maximum value of residual stress increased from 7.446 × 10 7 N/m 2 to 8.480 × 10 7 N/m 2 during the duration of a single short-circuit shock, and they all appeared on the disks of the primary winding near the end.The distribution of residual stresses is similar to the distribution of magnetic field leakage under short-circuit conditions as well as the distribution of short-circuit EF.
In the radial direction, the residual stress is concentrated between the two windings near the gap.In the axial direction, it is concentrated in the middle of the windings.The primary winding is subjected to a larger radial short-circuit EF than the secondary winding, so the overall residual stress value is higher than that of the residual stress of the secondary winding.In addition, the minimum value of residual stresses exists and is not zero because of the restraining effect of preload and gravity on the winding in the steady state.

Cumulative Effect of Plastic Strain under Multiple Consecutive Short-Circuit Shocks
Due to the multiple reclosures of power systems, recurring short-circuit shocks are a realistic risk to transformers.In this subsection, the stress-strain results obtained from a single short-circuit shock are used as the initial values for another short-circuit shock, and the cumulative effect of the winding under five successive short-circuit shocks is calculated.The distribution of winding plastic strain and residual stress after multiple successive short-circuit shocks is shown in Figure 15, and the results at t = 0.01 s after the fifth short-circuit shock are used as an example for comparative analysis.
Comparing the results in Figure 14a with those in Figure 16a, the plastic strain maximum increases by 1.601% after several short-circuit shocks, which is not negligible.Similarly, comparing the results in Figure 15a with those in Figure 16b, the maximum value of residual stress in the winding at t = 0.01 s after five short-circuit shocks increases to 1.088 × 10 8 N/m 2 .The minimum value also increases, which indicates that the initial residual stress value obtained from the steady-state calculation is related not only to the initial state force, such as the preload force, but also has a certain cumulative effect in the presence of successive short-circuit shocks.The distribution of residual stress increases with the number of short-circuit shocks and is influenced by the distribution of short-circuit EFs, and it is approximately the same as that of the first short-circuit shock.
lated.The distribution of winding plastic strain and residual stress after multi sive short-circuit shocks is shown in Figure 15, and the results at t = 0.01 s aft short-circuit shock are used as an example for comparative analysis.
Comparing the results in Figure 14a with those in Figure 16a, the plastic s imum increases by 1.601% after several short-circuit shocks, which is not negli ilarly, comparing the results in Figure 15a with those in Figure 16b, the maxim of residual stress in the winding at t = 0.01 s after five short-circuit shocks in 1.088 × 10 8 N/m 2 .The minimum value also increases, which indicates that the in ual stress value obtained from the steady-state calculation is related not only to state force, such as the preload force, but also has a certain cumulative effect i ence of successive short-circuit shocks.The distribution of residual stress incr the number of short-circuit shocks and is by the distribution of sh EFs, and it is approximately the same as that of the first short-circuit shock.
The results of plastic strain and residual stress of winding at different n short-circuit shocks are presented in Figure 17.The results of plastic strain and residual stress of winding at different numbers of short-circuit shocks are presented in Figure 17.
single short-circuit shock are used as the initial values for another short-circuit shock, an the cumulative effect of the winding under five successive short-circuit shocks is calcu lated.The distribution of winding plastic strain and residual stress after multiple succe sive short-circuit shocks is shown in Figure 15, and the results at t = 0.01 s after the fift short-circuit shock are used as an example for comparative analysis.
Comparing the results in Figure 14a with those in Figure 16a, the plastic strain ma imum increases by 1.601% after several short-circuit shocks, which is not negligible.Sim ilarly, comparing the results in Figure 15a with those in Figure 16b, the maximum valu of residual stress in the winding at t = 0.01 s after five short-circuit shocks increases t 1.088 × 10 8 N/m 2 .The minimum value also increases, which indicates that the initial resid ual stress value obtained from the steady-state calculation is related not only to the initi state force, such as the preload force, but also has a certain cumulative effect in the pre ence of successive short-circuit shocks.The distribution of residual stress increases wit the number of short-circuit shocks and is influenced by the distribution of short-circu EFs, and it is approximately the same as that of the first short-circuit shock.
The results of plastic strain and residual stress of winding at different numbers o short-circuit shocks are presented in Figure 17.According to the trend of winding plastic strain and residual stress with an increasing number of short-circuit shocks in Figure 17, the residual stress and plastic strain of the winding rises with the increase in the number of short-circuit shocks.The maximum increase in plastic strain and residual stress arises from the second short-circuit shock, and then it accumulates slowly with the increase in the number of shocks.

Conclusions
The vibration and deformation of transformer windings during short-circuit shocks are generally calibrated statically with the maximum value of the short-circuit current.Compared with most of the existing literature using 2D axisymmetric analysis, in this work, we carry out the dynamic analysis of a 3D model, in which the time evolutions of the short-circuit current and short-circuit electromagnetic force are taken into account.Transient 3D analysis allows flexible geometry and is able to incorporate the nonlinear material properties easily.
Based on our analysis regarding a 110 kV power transformer, it is found that the axial EF tends to compress the winding towards the middle.At the two ends of the windings, the amplitude of axial EF can reach its maximum of 4 × 10 7 N, while this value decreases to nearly zero for the middle disks.This is understandable because the radial leakage exhibits a similar distribution, namely, the peak amplitude appears at the two ends.
The axial vibration displacements of the winding under short-circuit conditions with varying pad area ratios and preload forces are investigated.It is concluded that increasing the area proportion of spacers in each layer, from 20% to 40%, can considerably improve the axial vibration stability, and a larger proportion is beneficial to improving the axial vibration stability.The vibration stability of LV and HV windings are also improved with the growth of the axial preload from 3 MPa to 6 MPa.In our example, a preload force of no less than 5 MPa can completely avoid instability.
The influence of aging temperature on the radial deformation of the windings is also analyzed.It is found that the maximum surface stress of the windings drops by 34% when the aging temperature varies from 20 • C to 130 • C, which indicates that excessively high aging temperature leads to potential axial instability accidents.
The nonlinear material properties of the insulating components and disks are taken into account, and the cumulative effects of multiple short-circuit current peaks under a single short-circuit shock as well as multiple consecutive short-circuit shocks are illustrated.It is shown that the plastic strain and residual stress dramatically increase by roughly 300% and 40.7% when the second short-circuit shock occurs, respectively, then they grow very slowly under the following shocks.

Figure 1 .
Figure 1.Equivalent model of transformers under short-circuit condition.

Figure 2 .
Figure 2. Finite element model of the transformer: (a) axisymmetric model and (b) 3D multilayer model.

.
The nonlinear elastic modulus is specified in Section

Figure 3 .
Figure 3. B-H curve of the core.
t = 0.01 s.The transient short-cir rent obtained by Equation (2) is shown in Figure4, and it is applied to the wind excitation source.

Figure 4 .
Figure 4. Transient short-circuit current in different windings.

Figure 5 .
Figure 5. Leakage flux of transformer windings under fault currents.

Figure 7 .
Figure 7. Different contact areas between disks and spacers.

Figure 7 .
Figure 7. Different contact areas between disks and spacers.

Figure 8 .
Figure 8. Difference in displacement of adjacent disks with numbers of (a) LV windings and windings.

Figure 8 .
Figure 8. Difference in displacement of adjacent disks with numbers of (a) LV windings and (b) HV windings.

Figure 9 .
Figure 9. Difference in displacement of adjacent disks with numbers of (a) LV winding and winding under 3 MPa, 4 MPa, 5 MPa, and 6 MPa preload.

Figure 9 .
Figure 9. Difference in displacement of adjacent disks with numbers of (a) LV winding and (b) HV winding under 3 MPa, 4 MPa, 5 MPa, and 6 MPa preload.

3. 3 . 3 .
Vibration and Deformation with Aging Insulated Spacers The support strength of spacers differs at different aging temperatures.Based on the axisymmetric model of the transformer, the corresponding Young's modulus expressions at different aging temperatures are imported as the material parameters of the model to simulate the axial vibration of the winding at 20 • C, 60 • C, and 130 • C. Comparing the displacement of the winding on each disk at different time instants, the distribution of the maximum value of the axial displacement difference of each disk is obtained, as shown in Figure 10.Energies 2023, 16, x FOR PEER REVIEW (a) (b)

Figure 10 .
Figure 10.Difference in displacement of adjacent disks with numbers of (a) LV winding and winding at 20 °C, 60 °C, and 130 °C.

Figure 10 .
Figure 10.Difference in displacement of adjacent disks with numbers of (a) LV winding and (b) HV winding at 20 • C, 60 • C, and 130 • C.

Figure 12 .
Figure 12.Stress-strain curve of the disk.

Figure 12 .
Figure 12.Stress-strain curve of the disk.

Figure 13 .
Figure 13.Stress-strain curve of the spacer.

Figure 14 .
Figure 14.Plastic strain distribution under single short-circuit impact at (a) t = 0.01 s and (b) t = 0.1 s.

Figure 15 .
Figure 15.Residual stress distribution under single short-circuit shock at (a) t = 0.01 s and (b) t = 0.1 s.

Figure 15 .
Figure 15.Residual stress distribution under single short-circuit shock at (a) t = 0.01 s and s.

Figure 15 .
Figure 15.Residual stress distribution under single short-circuit shock at (a) t = 0.01 s and (b) t = 0.1 s.

Figure 16 .
Figure 16.The distribution of winding (a) plastic strain and (b) residual stress after mul sive short-circuit shocks.

Figure 17 .
Figure 17.The distribution of winding: plastic strain and residual stress for various short-circuit shocks.

Figure 16 .
Figure 16.The distribution of winding (a) plastic strain and (b) residual stress after multiple successive short-circuit shocks.

Figure 16 .
Figure 16.The distribution of winding (a) plastic strain and (b) residual stress after multiple succe sive short-circuit shocks.

Figure 17 .
Figure 17.The distribution of winding: plastic strain and residual stress for various numbers short-circuit shocks.

Figure 17 .
Figure 17.The distribution of winding: plastic strain and residual stress for various numbers of short-circuit shocks.