Application of Duality-Based Equivalent Circuits for Modeling Multilimb Transformers Using Alternative Input Parameters

The principle of duality is applied for electromagnetic transient (EMT) modeling of industry scale (i.e. 50, 390 MVA) multilimb transformers. While saturation, hysteresis, deep-saturation, and remanent flux are accounted for, the need for transformer internal design information such as core dimension or material is eliminated. This is achieved by formulating the equivalent circuits with an alternative set of parameters that are either provided by the manufacturer or can be determined using conventional techniques. Open-circuit tests confirm that the models produce accurate excitation currents at different saturation levels when compared with measurement results. Furthermore, the models facilitate correct short-circuit condition with support for arbitrary number of windings. Upon validating the models, inrush current is simulated and the worst-case scenario is determined due to potential remanent flux values. The findings agree with an established EMT simulation model as well as manufacturer analytical approximations. Simulated hysteresis loops are also investigated.


I. INTRODUCTION
In recent years, the importance of sustainable development has urged the electric power generation industry to rapidly expand investments in environmental-friendly forms of generating electric power such as solar and wind power technologies. Therefore, interest have grown in using power system equipment that facilitates generation, transmission and distribution of electric power resulted from renewable energy. This includes power electronics, solar photovoltaics, inverter transformers, etc. Moreover, grids which carry renewable energy typically consist of a combination of small-and medium-sized inverter-based generators such as solar panels and wind turbines as well as energy storage systems such as batteries and flywheels. As a result, the use of small-and medium-sized power system equipment has The associate editor coordinating the review of this manuscript and approving it for publication was Mehmet Alper Uslu. generally increased in recent years. One important example of such equipment are three-phase transformers which are constructed on a single multilimb core as it is economically and logistically more advantageous to build small-and medium-sized transformers on a single core. Such cores are typically 3-legged or 5-legged as depicted in Fig. 1. Consequently in recent years, there has been an increased demand for more accurate and numerically robust three-and five-limb transformer models for wide-band [1] and low-frequency [2] electromagnetic transient (EMT) analysis. Among the available techniques for the latter, the principle of duality has been shown to be accurate and numerically robust for modeling many transformer topologies including multilimb cores [3]- [9]. Such models require internal design information (e.g. core material, core dimensions, number of turns in the windings) or results from unconventional tests (i.e. with open terminals) in order to determine parameters of the magnetizing branches. These parameters are pertinent to open-circuit test condition and heavily influence saturation studies. This is a practical limitation as such information is seldom available. Furthermore, to the best of our knowledge, available multilimb transformer models are limited to 2-and 3-winding cases. This is perhaps due to the difficulties arising in the enforcement of zero-sequence leakage inductances for arbitrary number of windings. Although in this paper we do not provide a solution to this problem and neglect zero-sequence tests, we attempt to clarify this difficulty for future studies. It should be noted that while transformers do not typically have more than three physical windings, inclusion of different tap positions with accessible terminals in EMT simulations requires transformer models to be extended beyond the physical number of windings and thus essential in some EMT studies.
In this paper, we formulate duality-derived models with an alternative set of required input parameters. They are namely; the nameplate data, factory acceptance test (FAT) report, core aspect-ratios, air-core inductance, hysteresis loop width, and the knee voltage. The resulting equivalent circuits fall into the category of gray-box models [3] as they are configured to match the target topology ( Fig. 1) with parameters to be determined. Since all the required parameters are traditionally used in EMT computations [10]- [13], they are either provided by the manufacturer or can be obtained via conventional calculation, measurement, or approximation techniques, some of which are reviewed in this paper (Section II-E). Consequently, the models exhibit accurate open-circuit test behavior. In order to ensure correct short-circuit test condition, the method of [14] is applied to multilimb core models with support for arbitrary number of windings. The resulting duality-based equivalent circuits are used to model a 50 MVA 3-limb 3-winding transformer and a 390 MVA 5-limb 2-winding unit. Results are compared against excitation current measurements under several saturation conditions as well as leakage inductance tests. To accommodate the study, sufficient saturation levels of the core have been achieved at the manufacturer's testing facility despite the industrial scales of the transformers. Subsequently, the models are used for calculating inrush current due to different remanent flux scenarios. Finally, the computed hysteresis loops are studied.

A. PRELIMINARIES
The first step towards forming the equivalent circuit of a transformer based on the principle of duality, is to determine the extent to which different parts of the transformer geometry are to be discretized. That is, how many magnetizing branches (linear or non-linear) are used to model different parts of the iron core and possibly tank, as well as how many linear elements are used to represent the leakage inductances. In this section, the principles that are followed in this paper are discussed. Throughout the paper, discussions are made for both 3-and 5-limb core models where winding limb (w) and yoke (y) apply to both cases, while outer limb (o) is only applicable to the 5-limb case.

1) IRON CORE
Different authors have proposed different approaches for modeling the iron core of multilimb transformers based on the principle of duality [3], [5]- [9]. Finer discretizations (i.e. using more non-linear inductors to represent the iron core) are believed to produce more accurate results. However, when forming duality-based transformer models in EMT-type programs, it is not only required to correctly represent the magnetic behavior of the device, it is also imperative to define the circuit parameters based on the available data. Fig. 1 depicts the approach taken in this paper for 3-and 5-limb cases. This is a generalization from linear inductors in [15] to non-linear inductors and linear resistors. That is, in accordance with the guidelines for slow transients [16], the non-linear (hysteresis) inductances model saturation, hysteresis, and deep-saturation phenomena, while linear resistances account for the core eddy current loss. Note that in Fig. 1 and throughout the paper, parameters with a tilde (e.g.L) have saturable/hysteresis nature, otherwise they are assumed to be linear (e.g. R). It was shown in [15] that discretizing the iron core of 3-and 5-limb transformers as shown in Fig. 1, allows for determining the linear magnetizing parameters. It does not require detail information about the core or unconventional tests. Instead, it ensures the accuracy of EMT simulations for the open-circuit condition based on the core aspect-ratios In (1), the cross-sectional area and the effective length of the winding limbs (s w , l w ), yokes (s y , l y ), and outer limbs (s o , l o ) are as shown in Fig. 1. In [15], it was shown that by knowing the magnetizing current I m along with r s and r l of a 3-limb core transformer, it is possible to compute linear inductances for winding limb L w and yoke L y with machine precision accuracy. The 5-limb core construction was also discussed where in addition to r s and r l , the core aspect-ratios for the outer limb (i.e. r s and r l ) along with I m could be used to precisely compute linear inductances for the winding limb L w , yoke L y , and outer-limb L o . In other words, by knowing the average magnetizing current and core aspect-ratios, one can exactly determine the associated linear inductances for individual core sections of 3-limb and 5-limb transformers. The same concept can be used in determining other linear core parameters from pertinent average values. Table 1 summarizes such parameters useful for the purpose of this paper. From Table 1, it is realized that one can obtain linear resistances R w , R y , R o from the average core eddy current losses I c . Also, the average air-core inductance L air can be used to determine air-core inductances for individual core sections L air w , L air y , L air o . Finally, the average hysteresis loop width D can be used to calculate loop widths associated with winding limb D w , yoke D y , and outer limb D o . Note that [15] contains the mathematical proof for the first row of Table 1 (i.e. computing L w , L y , L o ) and the other three rows can be proven similarly. As shown in Section II-C3, these parameters can be used to formulate non-linear magnetizing branches and thus eliminating the need for internal design information or unconventional tests.

2) LEAKAGE INDUCTANCES
Several methodologies exist when dealing with the leakage inductances of duality models [2]. In this paper, we adopt the mutual couplings approach originally introduced in [8] for 3-winding transformers and generalized to multiwinding transformers in [14]. It is simple and general while being accurate and numerically robust. In addition to the already published work [2], [8], [14], the authors have verified it to be reliable for leakage inductance modeling where up to 12 windings are assumed. Therefore, it is suitable for EMT studies involving multiple tap positions. This includes 3-phase bank as well as 3-limb and 5-limb transformers studied in this paper.
Despite the aforementioned advantages, it should be noted that models based on the method of [14] (including the ones presented in this paper), have two limitations. First, they are inherently ''non-reversible'' as studied in [18]. Reversible and non-reversible models behave similarly in the linear region and around the knee area. However, in the deep-saturation region, non-reversible models are less accommodating. That is, they are tuned to produce correct deep-saturation results (e.g. inrush current) when energized from a specific winding. Parameter adjustments may be required if energized from different windings. The adjustments are primarily needed for the air-core inductance as deep-saturation results are mostly dependent on this parameter. In contrast, reversible models produce accurate results regardless of which terminal they are energized from [6], [17], [18]. This is achieved by incorporating air-core inductances seen from all terminals into the formulation which typically results in a system of nonlinear algebraic equations. Therefore, reversible models are more computationally expensive and require extra input parameters compared with their non-reversible counterparts. Considering the fact that most transformers are meant to be energized from a unique terminal, non-reversible models are commonly preferred unless power flow needs to switch direction during operation. Examples include transformers in a storage system and high-frequency transformers in dc-dc converters of dc systems. 1 The second limitation of using the method of [14] is that while the resulting n-winding models can enforce n(n − 1)/2 positive-sequence leakage inductances, they do not necessarily ensure all n(n − 1)/2 zero-sequence leakage inductance measurements. For 5-limb core as well as three-phase bank models, this limitation is insignificant as the zero-sequence flux closes its path mainly through the iron core. However, for 3-limb models, considering zero-sequence tests may become important, especially in studies involving unbalanced energizations or unbalanced faults [3], [5]. A 3-limb transformer model that ensures both positive-and zero-sequence leakage inductances with support for arbitrary number of windings will be a subject of future work.

B. EQUIVALENT CIRCUITS
The equivalent circuits shown in Fig. 2 are applied for modeling multiwinding 3-and 5-limb transformers. The mutual inductances for each phase enforce positive-sequence leakage measurements and allow for arbitrary number of windings i = 1, . . . , n. As typically done in duality-derived models, ideal transformers are used to isolate the core from the winding terminals where desired connections (i.e. delta, Y, Zigzag) can be configured. The turns ratios for ideal transformers (N i :N c ) are dictated by the single-phase winding voltage V i and the reference voltage used at the core (i.e. core voltage V c ). In Section II-D, guidelines for choosing V c are discussed. The parallel hysteresis inductor-resistor pairs in Fig. 2, are used to represent different parts of the iron cores consistent with the discretizations in Fig. 1. In both equivalent circuits, it is noticed that R w ||L w pairs are connected to the mutual inductances where winding 1 is connected. Other windings are separated from the iron core through mutual inductances. Therefore, as indicated in the figure, the specified air-core inductance L air is seen from winding 1 only. This makes the models non-reversible as discussed earlier in Section II-A2.
It is important to note that while the equivalent circuits of Fig. 2 are somewhat similar to existing low-frequency multilimb models such as [3], [4], [7]- [9], they have been adequately simplified according to the principles discussed in Section II-A. Such simplifications facilitate formulation with the available set of input parameters while providing support for arbitrary number of windings.

C. FORMULATION 1) LEAKAGE INDUCTANCES
The leakage inductances are typically measured by performing the standard impedance voltage tests [19]. It is based on short-circuit test between windings i and j, while all other windings are open. Therefore for an n-winding transformer, n(n − 1)/2 short-circuit tests are carried out. The supplied voltage V i at winding i is increased until about 1 per-unit (pu) current flows through winding i (i.e. I i 1 pu). As winding j is short-circuited, a small voltage at winding i is expected to produce 1 pu current (i.e. V i 1 pu). Subsequently, by neglecting winding resistances, the leakage inductances are computed as where ω = 2πf with f being the transformer's operating frequency. The leakage inductances L i,j in (2) are typically provided by the manufacturer via standard FAT report or nameplate as leakage reactances X i,j in pu. Note that in terms of real values, inductance L and reactance X are related as a function of frequency X = jωL. However, when pu values are considered, they are numerically equivalent and interchangeable. From a total of n(n − 1)/2 leakage inductances, n − 1 self inductances L i,i+1 can directly be used in the equivalent circuits of Fig. 2. The mutual inductances M i,j are obtained as [14] The generality of the above described leakage inductance modeling allows for splitting the physical windings to multiple sub-windings with appropriate voltage levels and leakage inductances. This is typically done when multiple tap positions are to be accessed during EMT simulations.

2) LINEAR RESISTANCES
The equivalent circuits of Fig. 2 consist of linear resistors. The copper losses for individual windings are represented using winding resistances (R 1 , . . . , R n ). The linear resistances R w , R y , and R o are included to model core eddy current losses in the winding limbs, yokes, and outer limbs, respectively. As summarized in Table 1, they can be determined using the method of [15] if the total core eddy current loss I c is known. As explained in Section II-E3, I c can be obtained from information in the FAT report or nameplate. Therefore, R w , R y , and R o can be computed.
It should be noted that incorporation of non-linear effects such as anomalous (excess) losses [20] is generally recommended, especially for high-frequency applications and could be introduced in this work too. Nevertheless, assuming linear core eddy current loss is known to be valid for many low-frequency transient applications [16]- [18] which corroborates with our experience including the results presented in this paper. Since inclusion of anomalous losses would require additional input parameters, in this work we have limited the core loss to linear core eddy current loss and hysteresis loss (see also Section II-E3).

3) HYSTERESIS INDUCTORS
The hysteresis inductors in the equivalent circuits shown in Fig. 2, represent the non-linear magnetizing currents corresponding to the different parts of the core. Using saturable or hysteresis inductors that would require non-linear characteristic of the core material and dimensions such as B-H data [6] or the Jiles-Atherton method [9] is applicable. However, in this paper our aim is to develop gray-box models that do not require such data. To do so, we start from the asymptotic equation used to model saturation of the core [13] (see also All parameters in (4) have already been defined except for the knee voltage k. In order to include the hysteresis effect,  (4) and (5).
the first line of (4) is extended to include the hysteresis loop width D [10] i q (λ) = ±( whereĩ 1 toĩ 4 are the hysteresis magnetizing currents for the first to fourth quadrants obtained from all four combinations of ± signs. Fig. 3 plots the saturation curve and major hysteresis loop defined in (4) and (5). Note that if the loop width is zero, (5) reduces to the asymptotic equations modeling only saturation with no hysteresis effect. Such representation of hysteresis/saturation is sometimes referred to as the ''basic hysteresis'' model [13]. Based on (5), the hysteresis inductance in all quadrants for winding limbL q w , yokeL q y , and outer limbL q o can be expressed as where the following definitions are used When obtaining the values of hysteresis inductances using (6) and (7), it is realized that L air α , L α , and D α , are available from Table 1. In other words, the air-core inductance (L air α ), equivalent linear inductance (L α at V = 1 pu), and the loop width (D α ) for individual core sections (α = w, y, o) are computed using the method of [15]. Subsequently, they are used in (6) and (7) for obtaining the hysteresis inductances associated with different parts of 3-limb and 5-limb cores. Moreover, except for λ α , other parameters (ω, V c , λ c ) are similar to that of (4) and (5) which can be computed from information available in the FAT report and/or nameplate. In order to derive λ α , we start from the classical definition of flux linkage based on the core and winding details [6] λ α = N · B · s α , α = w, y, o (8) where N is the number of turns of the energized winding and B is the magnetic flux density. In (8), s α refers to the core cross-sectional areas for winding limb s w , yoke s y , and outer limb s o as depicted in Fig. 1. Despite simplicity of (8), in many studies, neither N , B, nor s α are available to system designers and thus (8) has limited application. Alternatively, similar to the definition of λ in the second line of (4), we can write λ w as where k w is the knee voltage corresponding to the winding limb. As explained later in Section II-E4, k w can be obtained or approximated from information provided by the manufacturer. Other parameters (V c , ω) are also known and thus λ w can be computed using (9). Subsequently, using (1) and (8), it is easy to establish the following relations This makes all parameters in (6) and (7) available. Thus, the sought-after hysteresis inductancesL q α can be evaluated. It is important to realize that computing the hysteresis inductances as formulated in (6), requires input parameters that are provided by the manufacturer or otherwise can be determined using feasible techniques as exemplified in Section II-E. Therefore, the practically-challenging task of measuring magnetizing parameters with open-terminal tests [7], [8] is circumvented while rarely-available knowledge of the core dimensions or material properties [6], [9] is not required.
Another important characteristic of the hysteresis inductances in (6) should be noted. They provide smooth non-linear functions as opposed to other forms of non-linear functions such as piecewise linear. In other words, the non-linear functions applied for formulating the magnetizing branches in Fig. 2 are smooth as the core saturates and transitions from linear, to saturation, and deep-saturation regions. This is consistent with the physical nature of the saturation and hysteresis phenomena. Evidently, by increasing the number of segments in piecewise linear methods, one can obtain functions that can accurately represent the non-linear nature of the physical phenomena. Therefore, hysteresis-aware techniques based on both types of functions can be shown to produce accurate results. However, this increases the computational complexity of piecewise-linear-based methods and can have negative impact on the simulation speed [23] and numerical stability [24]. Therefore, smooth non-linear functions (such as the ones presented) are preferable in EMT simulations.

D. IMPLEMENTATION
One of the main advantages of duality models is the fact that they can be effortlessly implemented using standard dragand-drop circuit elements available in EMT-type programs. This is recognized by other authors [2], [6], [8], [14]. Additionally, it is important to realize that electric circuits formed from standard circuit elements available in EMT-type programs are likely to be robust against numerical inaccuracies and instabilities. This is because standard circuit elements are designed for efficient and reliable EMT simulations. Therefore, electric circuits based on such elements are preferred not merely for their ease of implementation but also due to their reliable numerical accuracy and stability under different operating conditions. In this work, the equivalent circuits of Fig. 2 are implemented in the commercial software PSCAD [13]. It provides circuit elements for all electrical components needed to form the equivalent circuits. When implementing duality models, an appropriate core voltage V c should be chosen. To avoid numerical instability, realistic turns ratios for ideal transformers are required and thus we use V c = Min(V i ) as suggested in [7]. The above implementation guidelines are followed in this work and the models are available in [25]. In Section III, it is shown through numerical results that building the equivalent circuits as suggested above, ensures accurate and numerically robust models. However, the recommended guidelines should not be considered as restrictions. As such, it is expected that other EMT-type programs are applicable and other reference voltages could be applied or explored.

E. PARAMETER DETERMINATION
The required input parameters used in formulating the presented equivalent circuits can be categorized into two groups. First, parameters that are readily available from standard FAT report and/or nameplate (e.g. winding voltages, leakage inductances, operating frequency, etc). Second, values which are not part of the standard FAT report or nameplate, but have traditionally been used in EMT studies. Therefore, they may be provided by the manufacturer (initially or upon request) or can be accurately estimated using already existing techniques. Using the first group is straightforward and requires no further discussion. The second group of parameters along with several useful estimation techniques are discussed next. The manufacturer may provide core dimensions in which case the ratios can be computed using (1). The ratios can also be approximated quite accurately from dimensionless schematics of the core. This is because only the ratio is needed rather than the actual dimensions. Furthermore, it may be possible to accurately approximate these ratios by merely looking at a dimensionless image of the tank using gray-box approximations such as [26]. A very simple but effective method for approximating 3-limb core aspect-ratios is worth noting; by knowing that 3-limb cores are designed to carry the same flux from winding limbs to yokes, s w is usually the same as s y and thus r s = 1. The other ratio can simply be approximated from a dimensionless 2-D image of the core or tank as r l = length/height.

2) AIR-CORE INDUCTANCE (L air )
The air-core inductance (like many other electrical parameters), may be accurately computed using rigorous computational electromagnetic techniques such as the finite element method [27] or the method-of-moments [28]. Analytical formulations are also available [29]. However, this requires detail information about the transformer inner design. Although manufacturers rarely provide such detail due to intellectual property concerns, they may provide L air if requested. Alternatively, measurement techniques such as [30] are available for accurate and safe measuring of the air-core inductance.
In some studies, rather than the actual air-core inductance of the transformer, a pre-determined inrush current is assumed for a certain condition and its corresponding L air is computed for further analysis. This is a common practice when performing sensitivity analysis specially if the transformer has not been built or purchased. Such ''reference inrush current'' may be due to field data from similar transformers, worst-case scenario calculations such as [31], [32], or even empirical values. 2 In such cases, it is possible to run multiple EMT simulations and experimentally adjust L air until matching inrush current is simulated. Since the inrush current is primarily impacted by the air-core inductance, other parameters can be left unchanged when running multiple EMT simulations. It is expected that the model behavior in the linear region and around the knee area would not be impacted by different L air values. Hence such experimental finding of the air-core inductance is a relatively easy process. Note however that inrush current is highly sensitive to the resistance of the source and the connection cable, the remanence, and the circuit breaker. Hence, care must be taken when setting up such simulation experiments. It is important to note that realistic air-core inductance values are typically in the range of 0.1 pu to 0.3 pu. Too much deviations from this range is a clear indication that L air has not been properly computed.

3) LOOP WIDTH (D) AND CORE EDDY CURRENT LOSS (I c )
Standard FAT reports include ''no load loss'' which is sometimes referred to as ''total core loss''. The dominating factors contributing to such loss near the fundamental frequency are the core eddy current loss I c and hysteresis loss I h . In the presented models, the core eddy current loss is considered using linear resistors parallel with hysteresis inductors. The hysteresis loss is dictated by the loop width D. Therefore, by knowing how much of the total core loss corresponds to the core eddy current loss, the ratio I c /(I c + I h ) is determined and can be used to set the value of the core eddy current loss in the model. Subsequently, D can be found by performing multiple EMT experiments until the reported no load loss is produced by the model. Similarly, if the loop width is provided by the manufacturer, one can set the value of D and find the corresponding I c via EMT experiments. In the case that the ratio I c /(I c + I h ) is not known, a 50% − 50% split between I c and I h may be assumed since other non-linear effects such as anomalous losses are neglected.

4) WINDING LIMB KNEE VOLTAGE (k w )
Unlike other parameters discussed in this paper, the knee voltage does not have a direct physical meaning. Rather, it is mathematically defined to form asymptotic equations. Nevertheless, the concept of knee voltage has traditionally been used for representing saturation/hysteresis of transformers' core [13]. Therefore, manufacturers may provide this parameter. For single-phase and three-phase bank transformers, only one knee voltage k is defined (4). For multilimb transformers, the provided knee voltage typically corresponds to winding limb k w , which can be used in (9). If k y or k o are available, they can be used to obtain k w using core aspect-ratios similar to (10). In the case that the knee voltage is not provided by the manufacturer, system engineers may rely on optimization techniques to approximate it from saturation characteristic curve (V -I data) [34]. Values in the range of 1 pu to 1.25 pu may be considered realistic but too much deviations from this range is a sign of numerical inaccuracy or inappropriate approximation mechanism.

III. NUMERICAL RESULTS
The transformers under consideration are designed and manufactured by a Winnipeg transformer plant, PTI Transformers LP and further commissioned in the relevant industry. Due to the industrial scales of the devices (i.e. 50, 390 MVA) and the nature of the study specially saturation experiments, tests have been performed at the company's testing facility [35] with additional safety measures.

A. TRANSFORMERS' DETAIL
Details of the 3-limb transformer are given in Table 2. The air-core reactance X air is associated with the 138 kV terminal. As the model is non-reversible, winding 1 is assigned to this terminal. The other two windings are arbitrarily assigned to 13.8 kV and 6.972 kV terminals. Note that the tertiary winding is buried and thus X 1,3 , X 2,3 as well as R 1,3 , R 2,3 are estimated values. The loop width is available as a percentage of the magnetizing current I m = I 2 − I 2 c . Tests have been carried out for the base power rating of 50 MVA.
A 2-winding 5-limb transformer is considered with details given in Table 3. The provided X air is for the 238 kV side. Thus, winding 1 is used for this terminal. Tests are performed for the base power rating of 390 MVA.
All parameters in Tables 2 and 3 are supplied by the manufacturer except for the knee voltage and loop width which are determined using the techniques explained in  Section II-E. While most of the information in the tables are needed to build duality-based models of Fig. 2, extra information is added for interested readers. This includes core dimensions (s w , l w , s y , l y ), nominal magnetic flux density B, number of turns N 1 , and core material (23ZDKH85 [36]) where details such as B-H data can be extracted.

B. MODEL VALIDATION 1) LEAKAGE INDUCTANCES
The three short-circuit tests for the 3-limb model are simulated and results are compared with the leakage reactances in Table 2. The simulated values are X 1,2 = 0.076055 pu, X 1,3 = 0.114039 pu, and X 2,3 = 0.136033 pu. In all cases, errors less than 0.1% is observed. For the the 5-limb case, the simulated leakage reactance is X 1,2 = 0.06843 pu which is accurate according to Table 3. Figs. 4-6 compare excitation current waveforms from simulation and measurement tests for the 3-limb case. The agreement between measurement tests and simulations may be considered satisfactory. This is despite the fact that some deviations between measurement and simulation results exist VOLUME 8, 2020  mainly due to measurement under stress. For example in Fig. 6, measurement results for later cycles in phase A are plagued with DC shift and the observed mismatch should not be considered as error in the simulations. Also in Fig. 6, larger discrepancy between simulation and measurement results at the peak values may be explained by the low resolution of the measurement recordings not capturing sharp spikes due to high saturation levels.

2) EXCITATION CURRENT WITH DIFFERENT SATURATION LEVELS
For the 5-limb case, the exciting current waveforms are not available. Therefore, we compare the average current magnitude values I = (|I a | + |I b | + |I c |)/3 as a percentage of the rated current. This comparison may be sufficient in many studies as typically information about the core saturation is only available in terms of a V -I curve. Fig. 7 plots the results. It can be seen that the simulated saturation curve closely follows the behavior of the measured excitation currents.
As demonstrated above for both 3-and 5-limb cases, the non-reversible nature of the models does not prevent them from producing accurate open-circuit test results when energized from a winding other than winding 1. This is consistent  with the discussion made in Section II-A2. Therefore, excitation current measurements can be used to validate the models regardless of the energized winding.

C. INRUSH CURRENT CALCULATION
Inrush current is simulated by energizing the transformer from the HV side (winding 1) with 1 pu voltage. Other windings are left unloaded. As discussed earlier in Section II-A2, the models are non-reversible and thus inrush current experiments with energizations from other windings require X air to be adjusted. To avoid current limitation by the system impedance, ideal source is applied. The breaker is closed at the zero-crossing of phase A which leads to the maximum current in this phase.

1) NO REMANENCE
The inrush current is assumed with no remanent flux in the core. In theory, it is possible to physically perform such test by demagnetizing the transformer before energization. However, such tests are usually not performed unless laboratory size transformers are considered [6], [9]. Reasons include safety concerns, required voltage level, equipment damage, etc. Therefore we resort to validation by numerical convergence. Fig. 8 compares EMT simulation results with time-steps 50, 10, and 1 µs for the 3-and 5-limb cases. It is clear from the figure that both 3-and 5-limb models have numerically converged to their unique solutions. This also demonstrates the numerical robustness of the models against   increased time-steps when inrush current is simulated. This prevents the models from becoming a time-step bottleneck which can significantly impact runtime in EMT simulation of large systems.

2) WORST-CASE REMANENCE
Different potential residual flux values are investigated according to [33]. The worst-case scenario is found to be (0.8,0,-0.8) pu for both 3-and 5-limb cases. Results are summarized in Table 4 in terms of the peak values. For comparison, the table also includes simulation results from the unified magnetic equivalent circuit (UMEC) model [12], [13] as well as the values provided by the manufacturer based on analytical calculations [31]. However, it should be noted that there  are differences in the modeling and calculation techniques presented in Table 4. For example, the analytical formulas [31] are developed based on single-phase transformers and empirical scaling factors are used to approximate worst-case inrush current. The UMEC model neglects the hysteresis effect and the presence of the buried tertiary winding in the 3-limb case. Moreover, the UMEC model requires manual incorporation of residual flux by injecting DC currents at the line level [37] whereas the presented duality-based models enforce the remanent flux in different sections of the iron core using (6). Nevertheless, there is a general agreement among the results and the findings may be considered satisfactory.

D. HYSTERESIS LOOP CALCULATION
The λ-i hysteresis loop for different core sections can be computed using (6). As an example, Figs. 9 and 10 plot the open-circuit λ-i results for the middle winding-limb of 3-limb and 5-limb transformers, respectively. Here, winding 1 is energized with V 1 = 1 pu in both 3-and 5-limb cases. For comparison, results from the Jiles-Atherton method are also included where the magnetizing branches in the equivalent circuits of Fig. 2 are formulated using the magnetic characteristic of the material [21]. While some deviations exist between the simulated hysteresis loops especially in the non-major loops and around the shoulder area, the overall behavior match to a reasonable degree for both 3-limb and 5-limb cases. This is consistent with the results shown in Fig. 11, where it is observed that the methods have produced matching excitation currents. Therefore, the disagreements in Figs. 9 and 10 may be considered tolerable.

IV. CONCLUSION
In this paper, application of duality-based equivalent circuits for EMT simulation of positive-sequence leakage inductance, excitation current, inrush current, and hysteresis loop pertinent to multilimb multiwinding transformers is presented. Saturation, hysteresis, deep-saturation, and remanent flux are modeled while the required input parameters are either provided by the manufacturer or can be reliably estimated. Several experiments confirm the accuracy of the models. This makes them attractive alternatives for existing low-frequency models particularly when only limited information about the transformer is available or arbitrary number of windings are required. It is important to note that the lack of considering zero-sequence tests in the 3-limb case may impose limitations in its application in unbalance operation, a subject that is left for future work.  Dr. Gole is a member of the Long Range Planning Committee of the IEEE Power and Energy Society. He is a Fellow of the Canadian Academy of Engineering. For his contributions to the modeling of Flexible Ac Transmission System (FACTS) devices, he received the IEEE Nari Hingorani FACTS Award in 2007. VOLUME 8, 2020