1. Introduction

The reliable operation of transformers is of great importance for a whole power system as part of the backbone grid. However, unexpected vibration of transformers in operation seems to be quite a threat to such considerations and DC bias turns out to be one of the most significant reasons for this.

Affected by DC bias, deformation of the core and the winding will conduce to anomalous vibration as a kind of fault or even damage to transformers both directly and indirectly which is quite difficult to be detected or diagnosed by convention electrical test [1-3]. Hence, it’s fairly required to conduct research involving mechanism or so of such vibration phenomenon of transformers and so far some certain achievements have been reported already by literature from both in and abroad [4, 5]. A kind of novel approach to identify the condition of cores the windings employing acoustic-vibration method was stated in [6] that the mechanical faults like screw looseness were detected via time-frequency domain analysis of acquired transient vibration signal from a transformer. Another method to identify the characteristic frequency of transformer vibration is stated in [7] that the vibration of the core and structure deformation of transformers was analyzed utilizing magnetostriction characteristics. Besides, the excitation current in case of DC bias was tested and analyzed as described in [8]. And the nonlinear distribution of eddy current as well as the power loss in three dimensions was computed in [9]. The excitation current and harmonic distribution as well as hysteresis loop measurement of both single and three phases transformers were studied as stated in [10, 11]. Moreover, the experiment on excitation current and circuit-magnetic field corresponding to DC bias was researched in literature.

However, efforts like mentioned are just start-up research on transformer vibration yet, especially that on the variation pattern of abnormal vibration with presence of DC bias. In this paper, the relationship between several parameters of DC bias and vibration characteristics of transformers was investigated by both simulation and experiments and as conclusion, the effects on the core and the winding from DC bias were derived. Specially, the influence on flux density distribution from DC bias was identified by software simulation, which was coupled with mechanics field further. Afterwards, the variation pattern of mechanical stress, tension as well as displacement distribution of the core and winding was analyzed combined with classical Newton mechanics principle. And finally, the validity of simulation model as mentioned above was verified by the test platform specially established [12].

2. Mechanism of Transformer Vibration

2.1 Mechanism of the core vibration

Basically, as component of the core, the silicon steel sheets tend to vibrate periodically with the excitation frequency as a result of magnetostrictive effect. [13] Note that the period of magnetostriction as mentioned is just half of the power frequency cycle and it means the fundamental frequency of such vibration comes to be 100Hz which is twice the power frequency [14]. Again it’s much more noteworthy that with the development of transformer manufacture process like the improvement of the core lamination or colligation of core limbs with yokes by latitudeless epoxy glass adhesive tape, the vibration caused by electromagnetic force in sheets joint or between sheet discs turns out to be much slighter than before [15]. And the vibration of the core can be mostly attributed to magnetostriction of silicon steel sheets in the case where the prestress of the core is sufficiently large and silicon steel sheets are bonded tightly enough. The nonlinear magnetostriction of the core and the length differences of flux path along inner and outer edge of the core result in higher harmonic except for fundamental frequency in vibration spectrum of the core. Besides, if those silicon steel sheets haven’t been laminated tightly or the core bends and deforms impacted by its own weight, magnetostriction and meanwhile the vibration can be much more remarkable [16].

2.2 Mechanism of the winding vibration

As actually current carrying conductors located in leakage magnetic field, transformer windings are to be acted on electrodynamic force originated from alternating magnetic field, which is proportional to the square of current carried and nearly all the spectral components concentrate on fundamental frequency point. It’s just such electrodynamic force as mentioned that contributes to mechanical winding vibration and this kind of vibration is also likely to be passed to the other components [17]. Then if either of the HV(High-Voltage) or LV(Low- Voltage) wingding compaction degrades after deformation, displacement or collapse, the height difference between the two winding can be enlarged gradually, which will aggravate the ampere-turn unbalance between windings and meanwhile electrodynamic force tends to increase affected by magnetic flux leakage. Consequently, the winding vibration will be intensified. As to the vibration signal, compared to that under normal cases, amplitudes corresponding to each frequency points are to be changed and this change can be more notable if the displacement, looseness and deformation of the winding are more severe [18].

2.3 Mechanism of transformer vibration in case of DC bias

While DC bias arises, the flux density excited by positive and negative half wave respectively of the applied voltage exhibits asymmetry that the core flux density increases tremendously during the half with the same polarity as that of DC bias while it decreases during the opposite half [19]. As magnetostriction rate tends to increase with flux density, the elongation of silicon steel sheets will enlarge at the same time, leading vibration amplitude of the core to be raised [20].

On the other hand, flux flowing inside the core will saturate within the half segment of excitation current with presence of DC bias, which reduces permeability of the core and increases leakage flux at the same time. Furthermore, excitation current itself will rise caused by such flux saturation and consequently winding current will increase This in turn increases electromagnetic force acting on the two adjacent sections of core seam as well as magnetostriction rate. Ultimately, transformer vibration arises and is intensified for such reasons as mentioned [21].

3. Simulation Model

Typically, transformer vibration in case of DC bias is supposed to be considered as a situation of multi-physical field coupling. As a kind of solution to it proposed in this paper, firstly electric field model was established by transient current source provided in electromagnetic field block in the multi field coupling computation software Comsol. Then the inductive current unit was selected and magnetic field model was built with its excitation source set as calculation results of the above electric field model. These two steps consists, the first turn of coupling process namely electromagnetic field coupling and its distribution was identified. Whereafter the structural stress field model was established by solid stress/strain section provided in structural mechanics block where the electromagnetic force data calculated from the first model was loaded. This came to be the second coupling process and by these two model computation the distribution of stress, tension and displacement was obtained eventually.

3.1 Mathmatical model of vibration source of the core and the winding

This mathematical model of magnetic-mechanical field coupling inside transformers can be described by the form of matrix as (1) considering the magnetostriction of silicon steel sheets and anisotropy of magnetization characteristic.

Where M and K are magneto elastic coefficient matrix and mechanical elastic coefficient matrix respectively; C and D are coupling matrixes of magnetic-mechanical interaction which satisfies the relationship as C=DT. Specifically, matrix D indicates the impact on magnetic field distribution from mechanical vibration while accordingly matrix C indicates the opposite impact; A is introduced magnetic vector potential; u is the vibration displacement; Je is the external current density; fV and fΓ are respectively forces on solid and boundary surface of transformer the core.

Components of each matrixes element in expression (1) are defined along orthogonal axes x, y and z. As mentioned the silicon steel sheets are anisotropic. Thus, while solving the equation with Newton-Raphson method, nonlinear processing is to be implemented in advance. However, such work might complicate simulation model calculation, calling for much more computation cost. Besides, flux density concentrates within x-y plane though, there is still flux distributes along direction z at lap gap of the core namely the position with intensive vibration. Generally, magnetic force along direction z is calculated supposing the poverty of the core to be linear and isotropic. And based on such simplified excitation model, the volume density of magnetic force on the core can be expressed as:

Where f is density of volume force; J is current density; B is magnetic induction intensity(namely flux density); H is magnetic field intensity; μ is dielectric permeability; τ is volume density of dielectric. The first item in expression (2) refers to Lorentz force; The second item refers to volume force of material; The third item describes magnetostriction

Transformer winding vibration is definitely forced vibration in case of DC bias and the value of this force is proportional to the square of current, together with the periodicity of electrical excitation, the electromagnetic force and displacement can be written as:

Where , and T is the period of excitation current; b is the proportional coefficient which electromagnetic force to current; n is the positive integer; Im is the amplitude of sinusoidal current; M is the winding mass; C is vibration damping coefficient; K is insulation stiffness coefficient; F is the electromagnetic force suffered by winding; z is relative displacement of winding; g is acceleration of gravity.

3.2 Electromagnetic field model

A 5kVA three-phase dry type transformer with both its HV and LV being 400V was selected as modeling object here in the paper. And it’s simplified as we focus on the investigation of connections between cores windings vibration with simulated DC bias. Thus the following electric field equation can be derived based on the simplified model namely half section of the entire transformer.

Where ε0 is dielectric constant in free space and it equals to 8. 85×10−12 F/M; εr is relative dielectric constant; σ is conductivity; Je is external current density; V is the potential which is also the dependent variable in this equation.

While solving the electromagnetic field, the physical model was coupled with equivalent circuit utilizing fieldcircuit coupling method The external circuit was selected in electric magnetic field settings, which is illustrated by Fig. 1. And then the calculated current density data is loaded as excitation by magnetic field block. The magnetic field equation in solution region is as follows.

Fig. 1.System simulation circuit diagram

Where, μ0 is permeability in free space and its value is 4π×10−7H/m; μr is relative permeability; σ is conductivity; Je is external current density; A is magnetic vector potential, also dependent variable of the above Eq. (5)

As mentioned, the field-circuit coupling method employed here makes it much more convenient to load DC source referring to DC bias, and several other factors such as winding resistor or so can be considered, which definitely makes simulation of real transformer vibration to be exact and faithful. Through such simulation model, the distribution of H and B can be identified under different circumstances of DC bias. Besides, the following constitutive relation described by (6) and (7) will be utilized while solving the model.

Where E is the electric field intensity.

3.3 Structure force field model

In an actual mechanical or structural system the mass and elasticity of materials are supposed to be treated as distributed parameters. Nevertheless, it’s quite difficult to solve a computation model with distributed parameters, So the structural force field model here was simplified as discrete system containing springs with lumped mass connected with dampers. In this way, the dynamic characteristics of these structures can be illustrated by N orders of matrix differential equations which can be solved conventionally by computer simulation.

When solving structural force field model by Comsol, the solid stress block was selected then and the basic equation in solving region is as expression (8) with variables (u, v, w).

Where m, k and ξ are mass matrix, stiffness matrix and damping coefficient matrix respectively.

The force condition and displacement field distribution can be identified by solving the above equation, where electromagnetic field is coupled with structural force field by substitution of B and J into it. Specifically, in Comsol the filler of simulation model was set as air since the simulation object is a dry-type transformer and the surrounding air press was set as standard atmospheric pressure. Besides, materials of the core and winding were selected as copper and silicon steel provided in material library in Comsol.

It’s noteworthy that in this paper the constitutive relation of silicon steel sheets was determined by the nonlinear H-B curve of a specified ferromagnetic material, which was different from conventional approach that is modifying weak solution form so as to simulate characteristics of magnetostriction. Specifically, the corresponding relation between H and calculated B from the first coupling process was identified by interpolation method.

4. Simulations Analysis and Experimental Study

4.1. Establishment of electromagnetic model

In this part, the magnetic field distribution of transformer without DC bias was obtained firstly. As comparison, distribution of that was calculated again in case of DC bias so as to figure out the impact on magnetic field from different extent of DC bias.

4.1.1 Magnetic field distribution when DC offset voltage UDC=0V

Fig. 2 illustrates the distribution of magnetic flux density (B) when DC offset voltage is zero. Subfigure (a) and (b) refer to calculated distribution at the instant 0.005s and 0.015s corresponding to positive and negative peaks of excitation current. The scale of arrows in Fig. 2 as below demonstrates the distribution of B directly. And when comparing Figs. 1 and 2, it turns out that the extreme value of flux density as well as magnetic field distribution during both the positive and negative half cycle of excitation current is equal actually.

Fig. 2.UDC=0V Transformer magnetic flux density distribution

Fig. 3.UDC=10V Transformer magnetic flux density distribution

4.1.2. Magnetic field distribution when UDC=10V

When DC offset voltage UDC=10V, distribution of B is unbalanced in its positive and negative cycle. While applied positive offset voltage in simulation, it is found that the maximum of flux density at 0.005s namely the instant of positive peak is higher than that at 0.015s, the instant of negative peak, also higher than the peak value without applying DC bias. Explained bymagnetostriction principle, the magnetostriction of the core with presence of DC bias as above will definitely increase.

4.1.3. Magnetic field distribution when UDC=20V

Fig. 4 demonstrates that as DC bias increases, the unbalance of flux density in positive and negative cycle will be much more notable, especially the increase of flux density extreme value. Consequently, the magnetostriction coefficient is higher and vibration will be more significant.

Fig. 4.UDC=20V Transformer magnetic flux density distribution

4.2. Analysis of structural force field simulation

4.2.1 Computation of stress. tension and displacement fields of the core

Transformer vibration when without load originates mostly from vibration of the core. So the core model was setup with open circuit of secondary side to simulate an unloaded transformer. Besides, as in a genuine transformer, the upper and lower yokes are fixed by fasteners, the corresponding upper and lower boundaries of the core in simulation model were fixed in the same way as well.

Figs. 5 and Fig. 7 are distribution of stress, tension and displacement field at the instants of 0.0025s and 0.0075s namely the two peak instants when UDC=0V. Similarly, Fig. 8 and Fig. 10 are the distribution of stress, tension and displacement field at the same instants when UDC=10V.

Fig. 5UDC=0V The core stress and deformation

Fig. 6.UDC=0V The core tension distribution

Fig. 7UDC=0V The core displacement field distribution

Fig. 8.UDC=10V The core stress and deformation

Fig. 9.UDC=10V The core tension distribution

Fig. 10UDC=10V The core displacement field distribution

So as another group of simulation where UDC=20V, the corresponding distribution of stress, tension and displacement calculated then are not demonstrated here. Researchers in Japan such as M. Mizokami have ever tested the vibration of uniformly distributed measurement points in each surface of a 3 phase(3limb) core under the flux density of 1.3T and 1.7T. And vibration characteristics are compared under different flux density. Then came their conclusion that vibration exists mostly in magnetic return path of core and the vibration direction is perpendicular to front surface of core; vibration of the front and rear surface is more dramatic than the other surfaces. By comparison between Figs. 5, Fig. 8 and test result given by M. Mizokami, they are totally in agreement that the max stress in core appears in magnetic return path, which proves validity of this simulation method.

4.2.2 Computation of stress, tension and displacement fields of windings

When installing transformer windings, enough prestress is usually applied to both the top and bottom ends of it to guarantee its stability. So the upper and lower boundary is confined in simulation. Since calculated result from simulation model is time-varying solution, all the calculation results demonstrated here was extracted from solution data corresponding to the instant at 0.01s. As the transformer vibration can be greatest in case of short circuit, in this section the stress, tension and displacement fields of windings were calculated with secondary side being short circuit.

Figs. 11 and Fig. 12 are distribution graph of stress, tension and displacement fields of primary and secondary windings of transformers when UDC=0V. And Figs. 13 and Fig. 14 are distribution graph of those parameters above when UDC=10V.

Fig. 11.UDC=0V Primary winding stress, strain, displacement field distribution

Fig. 12.UDC=0V Secondary winding stress, strain, displacement field distribution

Fig. 13.UDC=10V Primary winding stress, strain, displacement field distribution

Fig. 14.UDC=10V Secondary winding stress, strain, displacement field distribution

So as another group of simulation where UDC=20V, the corresponding distribution of stress, tension and displacement field calculated then are not showed here. By comparison the displacement field distribution of primary and secondary winding under different DC bias, it’s found that with increase of applied DC offset voltage, the peak value of displacement field also increases significantly. For instance, as we can see from the simulation result, when UDC is increased by 20V, the amplitude of primary winding displacement is increased by 17% while that of secondary winding does not change significantly as it’s in leakage magnetic field with little flux density.

4.3 Experimental study

So as to verify validity of simulation model as mentioned above, a test platform of transformer vibration was built in this paper. It is consisted of a front-end acceleration sensor LC0166T, a signal conditioner LC0201-3, a data acquisition card YANHUA PCI-1714UL, man-machine interface based on virtual instrument and background databank monitoring system.

The experimental transformer was selected as a threephase dry type transformer. Its nominal capacity is 5kVA; voltage level is 400/400V; no-load current I0 is 12%; noload loss and short circuit loss are 75W and 160W; impedance voltage is 3.5%. In main magnetic circuit, the segment near silicon steel sheets vibrates much more significantly and this vibration of windings on iron beam can be measured exactly at such positions. So in experiments, three acceleration sensors were mounted on the outer surface of iron beam corresponding to phase A, B and C as illustrated by Fig. 15. In order to verify simulation and experiments mutually, the comparison and analysis as follows were implemented in case of UDC=0V.

Fig. 15.The test platform of transformer vibration

4.3.1 Validation of core vibration simulation model

In Fig. 15, vibration of transformer originates mostly from vibration of the core when without load. Thus the vibration acceleration curve of the point in magnetic return path with the greatest stress of a no-load transformer were extracted from structural force field simulation. The vibration waveform and spectrum are measured by acceleration sensor as demonstrated Also the corresponding spectrogram are plotted as illustrated by Fig. 16.

When compared the spectrograms in (a), (b), (c) of Fig. 16, their peak values appear at the frequency points of 200Hz, 300Hz, 400Hz and 500Hz in both simulation and measured results. The difference is that those spectrograms from simulation contain much higher harmonic. This is because the acceleration curve from simulation contains much burr, accordingly results in lots of high frequency components. However, we usually focus on the low frequency part and the frequency band over 1kHz tends to be neglected when analyzing transformer vibration. So generally, the simulation model proposed in this paper is proved to be of validation to simulate situations of transformer vibration definitely.

Fig. 16.The compare between core simulation model and measured by acceleration sensors

Regarding to the burr in simulation as mentioned above, it’s mainly caused by the higher frequency components in acceleration signal which arise in the coupling between magnetic field and mechanics field. The reasons are listed as follows.

(i) The magnetostriction of cores is nonlinear. (ii) The windings in transformers are modeled as coils, which causes differences between stiffness coefficients of two adjacent winding discs and two ends. And the stiffness coefficients are nonlinear consequently. (iii) The simulation step leads to burr in acceleration curves or so called the data mutation. When decreasingthe step size, the simulated curves will be smoother with fewer burrs. However, shorter step size calls for more computer memory at the same time which cannot be employed actually.

On the other hand, it is found that measured vibration signal also contains higher frequency components except for the interesting 100Hz component. Again, validation of simulation model is verified.

It’s noteworthy that 100Hz dominates the entire spectrumin simulation model meanwhile this dominance is not significant in measured data. This is because thecore is modeled ideally in this paper, which is free fromthe near mechanical components like fasteners, just rendering spectrum distribution in agreement with what mentioned in section 2.1 that spectral components concentrates near 100 Hz. While actually in a transformer, vibration data is measured with mechanical constraints on the core and resonance together with other mechanical components, which leads to the absence of 100Hz dominance. The model proposed just reflects essential mechanism of transformer vibration and is in agreement with actual case pretty well.

4.3.2Validation of winding vibration simulation model

As mentioned before, in case of short circuit, vibration of transformer is resulted from winding vibration. So by comparison between acceleration data solved in structural force field simulation and measured data, validation of winding vibration model can be verified. The vibration acceleration curves of windings and spectrograms are showed as Fig. 17 which measured by accelerometer B turns out to be much more significant than that from accelerometer A and C as demonstrated by Fig. 17 when setting secondary side of transformer be short circuit in simulation.

Fig. 17.The compare between core simulation model and measured by acceleration sensors

It’s found that simulated curve contains much burr as well, while it’s similar to measured curve as a whole. And Burrs with comparative high amplitude appear at the position of simulated curves corresponding to the bulges on measured curves. Frequency band of vibration acceleration ranges from 100Hz to 500Hz and frequency distribution concentrates on the band of 100Hz, 200Hz and 300Hz with mere higher frequency. These higher frequency components originate from burr on simulated curves as discussed already in simulation of the core vibration in last section. Again, simulation is proved to be consentaneous with measured data and of validation generally.

5. Conclusions

In this paper, the mathematical model of cores and windings vibration excited by electromagnetic force in transformer with presence of DC bias was proposed on the basis of their vibration mechanism. The model proposed was applied to the modeling of coupled field in transformers. Besides, the mathematical model of flux density distribution inside transformers was built in which the variation pattern of inner flux density during positive and negative cycles of excitation with DC bias is investigated. It’s found that distribution of flux density turns out to be quite different during the two polarity cycle. Whereafter, the structural force field with considerations of magnetostriction characteristics of ferromagnetic material was coupled with the electromagnetic field founded already, which was employed to study the impact on stress, tension and displacement distribution from DC bias. Furthermore, vibration acceleration curves of transformer the core and windings were plotted, which proved the validation of simulation model by comparison with the measured data from experiment platform ultimately.