Further Improvements in Topological Transformer Model Covering Core Saturation

This paper is devoted to improvements of the low-frequency topological model of a three-legged stacked-core transformer. It is shown that B-H hysteresis loops employed in the transformer model and ψ– i curves measured at transformer terminals are quite different in shape. A distinction is made between transformer modeling at moderate and deep saturations, separated conditionally by the level of technical saturation (near 2 Tesla) typical for grain-oriented steels. In the former case, in addition to the magnetic coupling of different phase windings through the core, an important part is played by air gaps at core joints. It is found that the effective gap length depends on the instantaneous magnetic flux and this dependence is proposed to implement by variable inductances. It is shown that regardless of the core saturation depth, an important role in the transformer modeling belongs to the distribution of the zero-sequence path between all the phases. In addition to accurate replications of the special purpose saturation test, the modeled results are in a close agreement with positive- and zero sequence data measured on a 50 kVA transformer, as well as with the measured inrush currents.


I. INTRODUCTION
The well-known problem of low-frequency transformer modeling is the reliable prediction of its behavior under high flux densities in the core. Different saturation levels of the legs and yokes necessitate the use of topological transformer models in which these core parts are represented by individual branches.
The difference in flux density levels of these parts is reproduced by shunting the wound legs of the model with individual reluctances which, taken together, constitute paths for the off-core (zero-sequence) flux flowing from yoke to yoke, partly directly and partly through the tank walls.
In contrast to inrush current events, when the legs are driven in the extra-deep saturation (often higher than 2.5 T), the relatively moderate core saturation under DC bias is characterized by flux densities not exceeding 2.0 T that is the value of the so-called "technical saturation" of grainoriented (transformer) steels [1], [2]. The moderate core saturation can be caused by geomagnetically induced currents or industrial impacts, such as DC-powered transportation systems.
Unlike inrush transients, transformer modeling in this saturation region has not received sufficient attention in the available literature. The problem is complicated by the impossibility of obtaining all necessary information from only measurements taken on the terminals of a three-phase transformer. This is because the waveform of each phase current depends on processes in all three wound legs [3]. To eliminate interdependencies of different phases, special measuring techniques were proposed in [3] where the rated voltage is only used to avoid the impeding influence of the transformer tank. The disadvantage of these techniques is that they do not distinguish processes in the legs and yokes that makes corresponding transformer models [3], [4], [5] not topological in the sense outlined above.
An appropriate solution was proposed in [6] where a special test was utilized to saturate the three-phase transformer by means of two increased single-phase voltages of opposite polarity [6]. The hysteresis loops measured at transformer terminals are dynamic (50-Hz) curves whose shape is determined in a complex manner by a number of interdependent factors. Among these are the influence of zero-sequence paths, the impact of core gaps and transformer capacitances, as well as the effect of core material properties. Eventually, the terminal loops do not relate to the legs and yokes taken separately. This complicates their direct use in a topological transformer model where processes in these core parts are considered individually. Therefore, the measured terminal loops are used in this work as reference patterns that should be reproduced by means of a specially fitted enhanced transformer model.
Because of the absence of reliable data about the core material and due to the always unknown transformer building factor (the ratio of the actual iron loss of a core to the estimated value based on the nominal loss of the core steel) [7], [8], a dynamic hysteresis model (DHM) with a set of predefined grain-oriented steels was used in the existing transformer model [9]. The need for its enhancement was caused by discrepancies between the regular shape of dynamic loops generated by the original DHM-based transformer model and the rounded hysteresis loops recorded in the measurements. The character of the rounding indicated the presence of variable air gaps at the core joints, while the negative slope of the measured cobratype loop near its waist pointed out the influence of transformer capacitances [10].
Details of the original experimental technique employed in the work are described in Section II, while peculiarities of the model are presented in Section III. Finally, using the model developed, different representations of the zerosequence path are discussed in Section IV.

II. EXPERIMENTAL SINGLE-PHASE SETUP
It is commonly recognized that the core material properties are of primary importance for modeling transformer transients [11]. If the core material and its hysteresis characteristic are unknown, the terminal measurements can be useful for obtaining relevant information. Most often, the flux linkage versus current hysteresis loops ( -i curves) are derived from the measured voltages and currents. However, due to the electromagnetic coupling of all phase windings through the core, the  -i characteristics derived from a three-phase no-load test are quite different [12], [9], and thus have a limited value for the subsequent modeling. For this reason, single-phase measurement techniques were developed in [3].
To mitigate the interphase coupling, an alternative single-phase measurement method was proposed in [6] where only inner, low-voltage (LV) windings are employed. In the setup depicted in Fig. 1, a three-phase transformer at no load is excited by two single-phase voltages (V test ) connected to the LV windings of the outer legs. Because of the 180° phase-shift between the voltages, the flux in the middle leg practically vanishes, thus isolating magnetically phase A and C of the core. This approach is tested on a three-legged two-winding transformer (50 kVA, 0.4/35 kV, YNyn0) depicted in Fig. 2 (transformer data are provided in the Appendix).  The LV transformer terminals A-X and C-Z are excited by two antiphase sinusoidal voltages formed by a 4quadrant power amplifier. These voltages are shown in Fig.  3(a), which also depicts the voltages measured across the LV terminals B-Y. Compared with the peak voltages of phases A and C (both are near 420 V), the maximum voltage of phase B (0.9 V) is almost indistinguishable from zero in the scale of Fig. 3(a). The latter justifies the assumption of the flux vanishing in the central leg.
To measure the  -i characteristic at an increased flux density level, the saturation test is carried out with perphase peak voltages of 420 V, which is some 30 percent higher than rated. Fig. 3(b) depicts the current waveforms, recorded during the saturation test. An asymmetry in the current waveforms in Fig. 3(b) is explained by high flux densities reached in the core. Even a small difference between voltages of phases A and C results in a noticeable difference between the peaks of corresponding currents.
To verify the model, the same test was also carried out for reduced peak voltages of 325 V.
Apart from the saturation test, winding capacitances to ground were measured on the transformer high-voltage (HV) side using the OMICRON Dirana test device. On an average, their values were about 3 nF. To avoid oscillations with model inductances, capacitances to ground were converted to the voltage of the inner windings and included on the LV side of the transformer model depicted in Fig. 4. Voltage and current waveforms recorded during the saturation test: a) Phase voltages across LV terminals; b) Current waveforms on the LV side.
The zero-sequence impedance was derived from the zero-sequence test according to [13]. At YNyn0 connection, all LV terminals of the unloaded transformer are connected together and supplied with a single-phase voltage V 0 . At the rms voltage V 0 = 52.73 V, the total measured current I 0 consumed by the transformer was equal to 70.5 A. For the subsequent modeling, the zero-sequence impedance was split up into three pairs of parallel-connected inductors (L 0A , L 0B , L 0C ) and resistors (R 0A , R 0B , R 0C ), as shown in Fig. 4. If these elements are the same for all three phases, then their values can be calculated using (1) and (2) through the threephase active (P 0 ) and reactive (Q 0 ) powers measured at frequency (f) of 50 Hz: The phase voltages and currents were measured with the DEWETRON DAQP-HV and the DAQP-LV together with Chauvin Arnoux MN38 current clamps.

A. MODEL STRUCTURE
In its basic features, the electric transformer model in Fig. 4 is similar to the ATPDraw model in [9]. The goal of this section is to improve this original model to make it capable of reproducing hysteresis loops measured in the saturation test.  The use of ATPDraw [14], which is a graphical, mousedriven preprocessor to the ATP version [15] of the electromagnetic transients program (EMTP), makes it possible to avoid the derivation of cumbersome equations and to reduce the study to the analysis of obtained illustrative graphs.
The legs and yokes are represented in the model by five DHM-elements, which are now included into the ATP library. The composite DHM consists of a historydependent static hysteresis model combined with dynamic components representing classical eddy-current and excess (anomalous) losses [16]. Starting with version 7 of 2019, the DHM can be found in the ATPDraw as L(i) Zirka-Moroz nonlinear branch [14].
Instead of the terminal  -i hysteresis loops, the DHM is based on B -H curves (dependences of magnetic induction on the magnetic field strength) taken from catalogs and extensive experiments conducted on the standardized Epstein frame [17]. At flux densities higher than 1.7 -1.8 T, ascending and descending branches formed by the model are merged into a single-valued saturation curve. At the level of 2.0-2.05 T (depending on the chosen material), the saturation B-H curve passes into a straight line with a slope equal to µ 0 = 4π•10 -7 H/m.
Due to the use of the DHM, the transformer model is able to replicate the measured no-load losses. The ratio of the actual iron losses to the losses obtained under standard conditions is characterized by the building factor, which usually ranges from 1.1 to 2.0 [7]. In the DHM employed, the loss control is carried out by varying coefficient K loss , which increases dynamic losses of the model.
The leakage representation includes the leakage inductances L 12 and inductances L 01 of the insulation channels between the core legs and inner LV windings. Because of the prescribed final slope of all saturation curves (dB/dH = µ 0 ) the value of L 01 is, in fact, the only fitting parameter of the model.
Three parallel connected L 0 and R 0 serve to portray the measured zero-sequence impedances. With the exception of Section IV, the values of these elements are considered to be the same in all three phases. Different L 0 -R 0 in Section IV are used to compare different representations of the zero-sequence path.
The 1:1 turns ratio of the ideal transformers (ITs) at LV terminals shows that the model parameters are referred to N 1 turns of the inner LV winding. The turns ratio n of the ITs at HV terminals is N 1 / N 2 where N 2 is number of turns in the HV winding. Winding resistances R 1 and R 2 are placed beyond the ITs, therefore their values remain unreferred. Measured transformer capacitances to ground were referred to the LV side and represented by capacitors C at the LV terminals.
In contrast to the majority of transformer models in which the air gaps in the core are neglected [10] or considered constant [9], the effective gaps at the core joints are found to be dependent on their instantaneous fluxes. Instead of evaluating gap characteristics using always uncertain geometry of butt-to-butt and interlamination air gaps, the resulting equivalent gaps in the transient model of Fig. 4 are modeled by means of nonlinear Type-98 inductors related to the corresponding legs. Neglecting flux fringing effects, the cross-section of the gap (S gap ) can be assumed equal to that of the leg (S leg ).
According to the saturation test, phases A and C of the no-load model in Fig. 4 are energized by equal antiphase voltages. For better model initialization and to reach symmetrical steady state, the voltages were increased from zero to their peak value of 420 V during 0.5 s. At the given cross-section of the core and the number of turns in the LV winding, peak flux densities in the legs A and C reach a quite high value of 1.931 T, while leg B remains unaffected.
The simulated current waveforms taken from the measuring switches SA and SC of the model are shown in Fig. 5 where the steady-state currents are close to the measured currents in Fig. 3(b).

B. MODEL IMPROVEMENTS
The main indicator of the model accuracy was the closeness of the calculated current peaks to the measured ones. With this in mind, some considerations can be useful before going to the detailed solution. The first is the drastic discrepancies between the measured magnetizing currents and those calculated with B-H curves provided by steel manufacturers [18]. Due to the impossibility of measuring the magnetizing current, it is better to discuss the input current (current i A in our case) and the hysteresis loop  A -i A , where the flux linkage  A is calculated conventionally by integrating the phase voltage V A (a small voltage drop over the winding resistance R 1 can be neglected to simplify the consideration).
It is useful initially to look at the terminal loop 1 in Fig. 6 calculated in the absence of core gaps (inductances L gap are infinitely large in this case and can be omitted in the model). By choosing a proper coefficient K loss of the DHM, namely K loss = 1.7, the width of this loop, as well as the width of all other calculated loops, is made close to that of the measured loop 2.
At the same time, loop 1 is more vertical in its upper part than the rounded loop 2. The tip of loop 1 is characterized by the current value of 18.5 A, which is markedly less than 24.4 A reached by the measured loop 2. The only mean to reach this current peak in the model with constant L gap is to decrease this inductance. As shown in Fig. 6, the loop with the desired current peak is obtained at L gap = 210 mH when the shape of the calculated loop 3 has nothing in common with the measured loop 2. Regarding the discrepancies arising from using manufacturer's B-H curves [18], it can be noted that quite different  A -i A loops 1 and 3 in Fig. 6 were obtained with the same B-H loop of the DHM inductors in the transformer model in Fig. 4. As an example, static and dynamic loops formed by the DHM inductor of leg A are depicted in Fig. 7. This shows that  -i loops measured at transformer terminals have only a vague similarity to the shape of B-H curves provided by manufacturers or catalogs. In particular, Fig. 6 illustrates a strong influence of the core gaps.
Although loop 3 in Fig. 6 is inclined too much to the right, its shape, together with the shape of loop 1, suggests the idea of the gap length Δ, which depends on the instantaneous flux in the core joint. The modeling above clearly shows that variable L gap should be large at small and moderate fluxes, while it has to be decreased when approaching saturation.
To explain this phenomenon, it can be considered that inductance L gap is inversely proportional to the gap reluctance R gap and hence to the gap length Δ [9]: Then, we can refer to the analysis in [18] where the small equivalent gap at small fluxes is explained by low reluctance paths of the magnetic flux around the large buttto-butt gaps. As the overlaps saturate, the flux is enforced to flow directly through the mentioned actual gaps causing a small L gap . Irregular sizes and positions of the core air gaps [19] made it impossible to estimate the gap length behavior directly. Besides, the negative slope of the measured loops near their waists (they can be seen in the inset of Fig. 8) pointed out the influence of winding capacitances. In this situation, it was found convenient to model the variable gap by the nonlinear inductance and fit its fluxcurrent curve  gap -i gap so as to approach the terminal loop measured at peak voltage of 420 V. The found curve  gapi gap in Fig. 9 was built by trials and errors method making use of ATP current-dependent inductor (Type-98) in which linear extrapolation of the last defined segment of the curve is used.
The loop calculated for the found dependence  gap -i gap and V A = 420 V is depicted by the solid line in Fig. 8. The accurate adjustment of its tip was achieved using L 01 = 0.7•L 12 . Here, the estimated factor 0.7 is close to its "standard value" of 0.5 [20] and to 0.62 used in [9]. The This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and The ratio  gap /i gap calculated for a given point of the curve is the current value of L gap . The value of Δ is then derived from (3) considering that S gap = S leg . As the air gap shares the same flux as the leg, the flux density B in the gap is calculated as  gap /(N 1 •S leg ). This also allows dependence Δ(B) to be built explicitly. In the inset of Fig. 9, it is represented using the logarithmic scale for axis Δ. It is remarkable that Δ(B) curve in Fig. 9 has qualitatively the same shape as the similar curve in [18, Fig. 4]. The simulated loops in Fig. 8 have resulted from opposite effects introduced by capacitance C and gap Δ that rotate the calculated loops counterclockwise and clockwise, respectively. It was found in numerical experiments that similar terminal loops can be obtained using different combinations of C and function Δ(B m ). As an example, function Δ(B m ) found for ten times smaller capacitance, i.e. for C = 0.3 nF, is shown in Fig. 10. The terminal loops calculated for these Δ(B m ) and C are similar to those shown in Fig. 8. It is therefore desirable to evaluate capacitance C first, and then adjust the corresponding flux-current dependence of the gap.
Together with the terminal loops in Fig. 6, the influence of the core gaps can be examined in Fig. 11 where waveform 2, calculated for the variable gap, is close to the measured current i a in Fig. 3(b). As pointed out above, curve 1 calculated for an ungapped core has underestimated peak, and its shape is different from the measured curve.  In concluding this subsection, we note the difference between the variable gap mechanism [18] implemented in this paper (it is inherent to laminated cores) and the constant average air gap technique leading to the effective permeability of a ferrite core [21], [22], [23].

C. MODEL TEST AT NO LOAD
The model parameters found by means of the saturation test were also used to predict transformer no-load currents and losses. Taking into account the asymmetry of the measured currents in Fig. 12(a), the overall coincidence of the measured and calculated currents in Fig. 12 is quite good. The same is related to a 9% underestimation in no-load losses that is well within 15% tolerance for component loss specified in the valid standard [13]. 7 VOLUME XX, 2022 Small discrepancies between the calculated and measured values are explained by experimental uncertainties and by the fact that "core gaps are different for each and every transformer limb" [3].

IV. COMPARISON OF DIFFERENT ZERO-SEQUENCE PATH REPRESENTATIONS
An important, but insufficiently studied issue in predicting transformer operations at high flux densities is a reliable representation of the transformer zero-sequence (ZS) path. As explained in Section II, it is formed by three individual L 0 -R 0 circuits shown in the model of Fig. 4. Such a distributed representation of the ZS path is exactly what makes the model topological -that is, one in which processes in the legs and yokes are described individually.
As an alternative, a simplified "concentrated" accounting for the ZS path is also widely encountered in the literature [3], [4], [5]. Due to the impossibility of measuring magnetization curves of legs and yokes separately, these core parts are merged into combined limbs A and C, that is, each of these lateral limbs consists of the leg and two adjoining yokes [3].
As there are no actual legs A and C in such a model, there are no ZS paths in parallel to these legs. This means that all ZS paths are collected together and placed in parallel to the central leg B. In terms of the model in Fig. 4, this requires that L 0B = 3•L 0 and R 0B = 3•R 0 where L 0 and R 0 are determined by (1) and (2). On the contrary, and to keep the structure of the model, its inductances L 0A and L 0C , as well as resistances R 0A and R 0C , have to be reduced to negligibly small values. As reduction factors of 0.01 and 0.001 give the same effect, inductances L 0A = L 0C = 0.001•L 0 and resistances R 0A = R 0C = 0.001•R 0 were used for definiteness in the model with the concentrated ZS path. With these parameters, the model was verified by simulating the standard ZS test, in which it showed the same behavior as the model with the distributed ZS path (the one with equal L 0 -R 0 in each phase).
The model with concentrated ZS path was used initially for simulating the saturation test. When the sources provide a peak voltage of 420 V, the model showed current peaks of 29 A that is some 20% higher than those recorded in the test and calculated with the model with distributed ZS path.
The observed result indicates that the model with the concentrated ZS path can provide significantly overestimated currents in lateral phases A and C during inrush events. This can be qualitatively explained using well-known expressions for inrush current peaks discussed in [20], [24]. The main component in the denominator of these expressions is the saturation inductance L sat , which is determined by the height h of the excited winding and its equivalent diameter d eqv : It was pointed out in [24] that the value of h in (4) does not include the length of the yokes, whereas in the model with the concentrated ZS path the yokes are included in the lateral limbs by definition. Such an artificial increase in h leads to a decrease in L sat and, accordingly, to an increase in current peaks.
Although (4) was derived for a single-phase transformer, the above explanation also applies to threephase units without delta-connected windings. To verify this conclusion, the unloaded transformer in steady-state was multiple times disconnected from the public grid and then reconnected to the grid. The goal was finding a case in which currents A and C would be maximum and approximately equal. The measured LV voltages in Fig.  13(a) and corresponding inrush currents in Fig. 13(b) were found suitable for the subsequent analysis.
By integrating the voltages in Fig. 13(a) and recalculating symmetrized flux linkages into flux densities, it was found that their values just before transformer reconnecting to the grid can be estimated as B   By entering these values into corresponding DHMinductors, the inrush currents shown by dashed lines in Fig. 13(b) are obtained for an appropriate initial phase of the grid voltage. The calculated and measured currents in Fig. 13(b) are almost indistinguishable.
As expected, the model with the concentrated ZS path behaves quite differently during inrush current events calculated for the same excitation scenario and residual flux densities. The current peaks in Fig. 14 are more than 30% higher than those in Fig. 13(b).
Thus, the inaccuracy of the transformer model with concentrated ZS flux path increases with the level of core saturation.

V. CONCLUDING REMARKS
This paper presents experience gained in modeling a threephase, three-legged transformer operating in a wide range of the core flux density. Particular attention is paid to using a topological transformer model under core saturation conditions. The correct topology of the model, that is its ability to distinguish processes in the legs and yokes, is ensured by distributing the zero-sequence (ZS) transformer impedance between three per-phase values.
The cost for using such a distributed ZS path model is the impossibility to measure magnetization curves of the legs and yokes separately. Because of this, as well as due to the usually unknown building factor, the transformer under test is modeled in ATPDraw [14] where legs and yokes are represented individually, making use of the dynamic hysteresis model (DHM) [16]. Among advantages of the DHM is the possibility to choose the core material from a set of predefined grain-oriented steels and the ability to reproduce arbitrarily high saturation.
To the best of the authors' knowledge, a distinction is made for the first time between transformer modeling at moderate and deep saturation levels. Accordingly, two different measurement techniques were utilized in addition to the standard tests to observe corresponding operations of a 50 kVA transformer available for detailed measurements. In the course of reproducing recorded terminal data, some changes to the existing transformer model became evident.
In a special purpose single-phase saturation test, the peak flux density in the core legs was found to be near 1.931 T. This is well in excess of 1.7-1.8 T typical for transformer cores assembled from grain-oriented steel, but markedly lower than 2.0 T representing the technical saturation. At this moderate saturation, a significant impact of core joints on terminal hysteresis loop was observed. It was found that the effective length of the core gap depends on the instantaneous magnetic flux in the adjoining leg. This finding led to the idea to represent variable air gaps by means of nonlinear inductances, which are easy to implement in the transformer transient model. This very feature has allowed the model developed to become a flexible tool for describing accurately the shape of the terminal hysteresis loop ( -i curve) whose tip determines the peak of the drawn current. VOLUME XX, 2022 The model fitted to the moderate current peak (24.5 A) measured in the saturation test was then verified by accurate prediction of inrush currents, which reached maxima of 586 and 692 A in the lateral phases A and C, respectively. The no-load currents and losses were also evaluated appropriately.
This allowed the distributed ZS path model to be used as reference in evaluating the acceptability of representing the measured ZS impedance by a single inductance 3L 0 and a resistance 3R 0 concentrated at the central core leg. The 30 % inrush current overestimation obtained using this "concentrated" model clearly illustrates its drawback at high flux densities. Speaking generally, this also shows that far not all of the models producing inrush currents are in fact reliable. An interesting finding is also that the impact of variable air gaps becomes negligible when modeling transformer in the extra-deep saturation. This is because the increase in the phase current due to the nonlinear gap (it is of the order of 20 A in the transformer considered) is negligibly small compared with that drawn by the transformer during the inrush current period. For example, the rounding of the terminal hysteresis loop seen in Fig. 8 is not perceptible in the loop in Fig. 15(a) calculated for the worst-case scenario regarding the current peak in phase A, Fig. 15(b).
It is worth recalling here that the total loss predicted by the DHM increases nonlinearly with the magnetization frequency and flux density [16]. The latter manifests itself in a wider area of the loop in Fig. 15(a) near saturation. This makes the calculated loop similar to the measured curve demonstrated in [25, Fig. 8].
All these observations and obtained results demonstrate that the accent in transformer modeling should be shifted from choosing specific core steel from the DHM library to details of the model itself. Among them is the loss coefficient K loss of the DHM, the air gap characteristics (at moderate core saturation), and the necessity of using distributed ZS path when the core is deeply saturated. The measurement-based approach presented in this paper provides a basis for reliable transformer modeling in a wide range of transformer operations.
Further work will focus on the model validation in a transformer back-to-back setup with superimposed direct current.
Appendix. Transformer data