Hindawi Publishing Corporation Mathematical Problems in Engineering Volume 2013, Article ID 215647,10 pages http://dx.doi.org/10.1155/2013/215647

Research Article

Simulations of Transformer Inrush Current by Using BDF-Based Numerical Methods

Amir TokiC,1 Ivo Uglesic,2 and Gorazd Stumberger3

1 Faculty of Electrical Engineering, Tuzla University, 75000 Tuzla, Bosnia and Herzegovina

2 Faculty of Electrical Engineering and Computing, Zagreb University, 11000 Zagreb, Croatia

3 Faculty of Electrical Engineering and Computer Science, Maribor University, 2000 Maribor, Slovenia

Correspondence should be addressed to Amir Tokic; amir.tokic@untz.ba Received 7 April 2013; Revised 8 July 2013; Accepted 9 July 2013 Academic Editor: Evangelos J. Sapountzakis

Copyright © 2013 Amir Tokic et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This paper describes three different ways of transformer modeling for inrush current simulations. The developed transformer models are not dependent on an integration step, thus they can be incorporated in a state-space form of stiff differential equation systems. The eigenvalue propagations during simulation time cause very stiff equation systems. The state-space equation systems are solved by using A- and L-stable numerical differentiation formulas (NDF2) method. This method suppresses spurious numerical oscillations in the transient simulations. The comparisons between measured and simulated inrush and steady-state transformer currents are done for all three of the proposed models. The realized nonlinear inductor, nonlinear resistor, and hysteresis model can be incorporated in the EMTP-type programs by using a combination of existing trapezoidal and proposed NDF2 methods.

1. Introduction

The transformer represents one of the essential elements in power systems. Proper modeling of the transformer is very important in different transient and steady-state power systems simulations. The nonlinearity of the transformer iron core is the most important parameter in simulations of low-frequency transients, such as transformer inrush current, ferroresonance, temporary overvoltages during transformer energizations, load rejections, and harmonic analysis. All these transients belong to a frequency range of up to 1 kHz [1] for which it is not necessary to take into account frequency dependencies of the transformer parameters. A standard Steinmetz T equivalent circuit model of a two-winding single-phase transformer is shown in Figure 1. In addition, there is an analogous n-shaped equivalent model. This model is rarely used due to the fact that the detailed transformer construction has to be known for determining the parameters of the n equivalent model.

Different types of T transformer models will be used in this paper in the first place. Transformer parameters can be divided into two basic groups: winding and iron-core parameters. In Figure 1, Rp and Rs represent the serial

winding resistances that generally include the Joule and eddy current losses in the windings, whilst Lp and Ls represent the serial leakage inductances divided amongst primary and secondary windings. The shunt branch parameters describe iron core behavior, where the resistance Rm represents the core eddy current losses, whilst inductance Lm represents the saturation or hysteresis phenomena in the transformer core. The winding parameters are linear, whilst the iron core parameters describe nonlinear phenomena during the analysis of low-frequency transformer transients.

In general, there are three different approaches when modeling transformer iron-core nonlinearities:

(a) linear resistor Rm and nonlinear inductor Lm,

(b) nonlinear resistor Rm and nonlinear inductor Lm,

(c) linear resistor Rm and nonlinear hysteresis inductor L m-

During the analysis of low-frequency transformer transients such as inrush current and ferroresonance, the more used model is the model (a) because of its simplicity; relatively rare used models are (b), and (c). For example,

t = T,

O n L-

NP : Ns Ls Rs

Figure 1: Standard T equivalent transformer circuit model.

>slj 's2 'sN

Figure 2: Nonlinear curve of core inductor.

[2-10] use the (a) based model of transformer, [11-14] use the (b) but [15-18] use the (c) based model of the transformer.

Typical procedures for obtaining an instantaneous nonlinear magnetizing curve and nonlinear curve of core losses, by the separation technique of losses from the magnetizing curve, are described in [19,20].

The modeling of nonlinear or hysteresis inductor and nonlinear core loss resistor can be done by using the curve fitting procedure for the approximation of nonlinearity or by using piecewise linear representation of the nonlinear curve.

The well-known curve fitting approximations are obtained by using the polynomial, exponential, arctg, and hyperbolic functions. The polynomial approximation was suggested by Mayergoyz in [21]. Exponential approximation was introduced by Trutt et al. in [22]. Arctg approximation was started by Karlqvist in [23], and the hyperbolic functions were introduced by Takacs in [24] and refined by Takacs in [25].

The piecewise linearization of a nonlinear curve for the modeling of nonlinear or hysteretic inductor, on the other hand, is used in [7-10]. Linearization of a nonlinear curve has certain advantages and disadvantages. Successful control of the numerical stability properties of the applied numerical methods is the main advantage [26-28]. Also, another advantage of this linearization is the improved speed compared to those classical methods used in calculations within nonlinear algebraic-differential equation systems. On the other hand, piecewise linearization has disadvantages related to overshooting effects [27].

2. Transformer Iron Core Modeling

The original way of modeling nonlinear transformer iron core in low-frequency transient analysis, using three mentioned models, is derived in this paper.

O_±_TYYYV

-TfYT\_

Lmk .O mk

Figure 3: (a) Nonlinear inductor and (b) equivalent model.

'rl 'r2 irM

Figure 4: Nonlinear curve of core resistor.

2.1. Modeling of Nonlinear Core Inductor. The nonlinear core inductor can be described using a magnetizing curve: magnetizing current versus magnetic flux, as shown in Figure 2.

When this curve is piecewise linear, the results are the input vectors magnetizing currents IL = [<'Si, iS2,..., iSN ]T

and magnetic fluxes Ol = [0 , 0S2,..., ]T, whereas N is the total number of magnetizing curve regions. The magnetizing current for the fcth linear region of the magnetizing curve, k = 1,2,..., (N - 1), is calculated using the equation

= + sign 0) ^ = + \ >

^ mk ^ mk

whereLmk = -$Sí)/(is+ -iSi) andISt = iSt-(1/Lm

sk+1 sk '

mk' Tsk'

Based on (1), it is possible to create a nonlinear inductor equivalent model using a linear inductor and a corresponding current source, Figure 3. This model is functionally dependent on the position of the operating point k. At the same time, contrary to EMTP-based equivalent models, the inductor model is not dependent on the integration step.

2.2. Modeling of Nonlinear Core Resistor. Analogous to the nonlinear core inductor, the nonlinear core resistor can be described using the curve: resistor current versus corresponding resistor voltage, as shown in Figure 4.

In this case, the input vectors are the resistor currents

IK., I

r2,... ,-rui and resistor voltages UR = [uri, T, where M is the total number of resistor curve

Rmk iRk

Rmk i r (b)

Figure 5: (a) Nonlinear resistor and (b) equivalent model.

regions. The resistor current for the fcth linear region of the resistor curve, k = 1,2,..., (M - 1), is calculated using the

equation

% = 7> Ur + sign (uR) Irk = ~n Ur + Sr,

where Rrnk = (urk+1 -Urk )/(irk+1 -irk ) and !rk = Kk-(1/Rmk )urk.

In the same way, it is possible to create a nonlinear resistor equivalent model using a linear resistor and a corresponding current source, Figure 5.

2.3. Modeling of Hysteresis Core Inductor. By definition, the major loop is the largest possible loop whose end points reach the technical saturation. Beyond the saturation points the hysteresis is a single-valued function. Symmetrical or asymmetrical minor hysteresis loops lie within the major loop. The following is assumed, based on experimental results, for hysteresis modeling (Figure 6) (a) the operating point trajectories lie within the major hysteresis loop, (b) in the direction of the flux increase (decrease), the operating point heads towards the lower (upper) branch of the major loop, and (c) after the flux reversal points are detected (1, 2, 3), the operating point trajectories will head towards the previous reversal point.

Therefore, following the symmetry, the major hysteresis loop is defined by the points of the lower (or upper) half of the major hysteresis loops, that is, by input vectors are: hysteresis currents ILh = [ih ,ih2,..., ihp] and magnetic fluxes <&Lh =

,..., <php ]T, where P is the total number of the major hysteresis loops.

The expression for the hysteresis current of the kth piece-wise region in terms of the main flux k = 1,2,..., (P-1), Figure 6, can be stipulated as

\ = J^^rn + sign {A$rn) Ihk > Lmk

where A$m = $m(tj) - <UVi), j = 1,2,3,.

($hk+1 - <1%)/(\+1 - \), and Ihk = \ - (1/L

In the general case, if it is assumed that the distance between the operating point and the major loop is a linear

Figure 6: Major and minor hysteresis loops.

function of the magnetic flux, for both the ascending and descending trajectory, the next equation is valid:

d = Vk$ + nk = sign (A^) 0 -

The unknown coefficients and qk are determined on condition that the two successive reversal points ($rev1, drev1) and ($rev2, drev2) lie on this line.

Based on relations (3)-(4) it can be proved that the expression for the hysteresis current of the kth linear region that describes the trajectory of the actual operating point is

+ sign (A&-(lhk -Mk), (5)

thk Lmt -AL

noting that ALmk = (sign(A$)pk/(sign(A$)pk - 1))Lmk and AIk = (1/Lmk )nk.

In this way the hysteresis inductor Lh in Figure 3(a) can be modeled by introducing a linear inductor Lh in parallel connection with a current source Shk, both functionally dependent on the actual linear region k, Figure 7 as follows:

L hk = L™k AL™k ' Sht = sign A)-(lht -AIk).

It is noted that all three developed models has a special routine for eliminating the possible overshooting effect [2628].

3. Inrush Current Modeling and Numerical Methods

The transformer inrush current is a well-known transient response to the energization of a transformer iron core. During transformer energization, depending on the value of the remanent flux, the magnetization curve and the breaker switching instant, the iron core magnetic flux can reach a twice higher value than the nominal operating value. In the case when the remanent flux is nearly as high as the saturation, then the inductance is reduced to near zero.

(a) Shk

-/YYYY_ Lhk '0 '

Figure 7: (a) Hysteresis inductor and (b) equivalent model.

The only initial impedance is the small ohmic resistance of the primary coil. As a result, the transformer iron core is driven into a saturation, and inrush current is produced. The transformer inrush current, which could be many times higher than the operational current, may cause irreparable damages to the circuit without taking into account the design stage.

The developed models for the nonlinear inductor, nonlinear resistor, and hysteresis inductor are very suitable for formulation state-space equations that describe low- frequency transformer transients such as transformer energization, Figure 8.

An arbitrary numerical method could be applied in such a state-space form. In the EMTP-type of programs, the elements are strongly dependent on the integration step; this fact becomes apparent when a trapezoidal numerical method is applied to the relevant branch. The general algorithm procedure for generating state-space matrixes is shown in [27]. According to this procedure, a state-space equation is developed that describes transformer transients during inrush current analysis:

dX(t) dt

= A jJcX(t) + BjJc U(t) = F(X,t),

where the state-space vector and input vector are:

X(t)=[iLp U(t)=[e(t) Sh Sr, ]T .

System matrixes are, respectively as follows:

A j,k =

r Rp + Rm¡

L PL mk

Bj,k =

— 0 L ,

— -Rm,

It is very important to discover whether the demonstrated system is a stiff system. Stiff systems are categorized as those different components of solutions which evolve on very different time scales occurring simultaneously, that is, the rates of change of the various components of the solutions which differ markedly. Stiffness is a property of system with strong implications for its practical solution using numerical methods. The stiffness of system is defined by two factors: stiffness ratio k and stiffness index k for all state matrixes A ;

max I Re

min I Re

= max J Re (A^

where i = 1,2,..., dim(Ajfc) are eigenvalues of all state matrixes Aj k.. In the case for k » 1 it isa stiff system, and in the case for k ^ >x> it is a very stiff system.

The state equations (7) are traditionally solved using the implicit A-stable, second-order trapezoidal method [7]:

Xn+1 = Xn + -[F (Xn, tn) + F (Xn+1 ,tn+1)],

n= 1,2,3,

The trapezoidal method is suitable for nonstiff and moderate stiff systems; however for a very stiff system this method could be inaccurate, and other techniques should be used [26-30]. Indeed, the trapezoidal method is not L-stable, such that, for state matrixes A^ with eigenvalues containing large negative real parts, this method produces unwanted numerical oscillations. In addition, the classical explicit numerical methods are not applicable for stiff or very stiff systems due to finite regions of numerical stability.

On the other hand, the Backward Differentiation Formulas (BDF_p) methods of _pth order, p < 5, as introduced by Gear in 1971, have the following form [29, 30]:

X -ymXn+i =AtF(Xn+i,tn+i).

The order of the truncation errors for the BDF_p methods is Atp+1 [29, 30]. The BDFp methods are A(a)~ and L-stable, with stability angles ranging between 90°, 90°, 86°, 73°, and 51° [25]. The L-stability properties of the BDF_p methods lead to the suppression of numerical oscillations,

Magnetizing current (A)

-*- Magnetizing curve —Hysteresis loop

Figure 9: Magnetizing curve and hysteresis loop.

whereas the trapezoidal method is not free from numerical oscillations [29, 30]. In addition, Numerical Differentiation Formulas (NDFp), as suggested by Klopfenstein in 1971, are a modification of the BDFp's methods with the next relations [31, 32]:

ti m (13)

+ Kyp (Xn+i -X^i),

where the coefficients are yp = Xm=1(1/m), the predicted

value is X<n+r1 = Xm=o VmXn, and parameter k was introduced by Shampine and Reichelt [32], which is a compromise made in balancing efficiency in step size and stability angle of the NDFp. These methods are also A(a)- and L-stable [32]. Compared with the BDFp's, the NDFp methods achieve the same accuracy as BDFp methods with a bigger step-size percents 26%, 26%, 26%, 12%, and 0%, respectively. This implies improvements in the efficiency of the methods. However, the percent changes in the stability angle are 0%, 0%, -7%, -10%, and 0%, respectively. It can be concluded the NDFp methods are more precise than the BDFp but not more stable. The other modifications of original BDFp

Resistor current (A)

- * - Linear resistor —b— Nonlinear resistor

Figure 10: Nonlinear core loss resistor curve.

Figure 11: Measured transformer inrush current.

methods are the extended BDF methods as proposed by Cash in [33], the generalized A-BDF methods as introduced by Fredebeul in [34], the diagonalizable extended BDF methods as suggested by Frank and van der Houwen in [35], a predictor modification to the extended BDF methods as introduced

H?: Power Analyzer MI7111 [Samples - Chart]

Samples

1,0 0,8 0,6 0,4 0,2 0,0 -0,2 -0,4 -0,6 -0,8 -1,0

Analysis Charts Settings

Steady state current (measurement)

: ............................

: / \ : / \

^y X Çy

x f i \ / : \ /

r i r i

11 —

14:37:47-Time

Figure 12: Measured transformer steady-state current.

Figure 13: Harmonic content of measured inrush current.

Inrush current (simulation)

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Time (s)

Figure 14: Inrush current simulation, Model 3, trapezoidal method, Af = 100 ^sec.

by Alberdi and Anza in [36], and the block BDF methods developed by Yatim et al. in [37].

The main aim of this paper was to explore the types of differential equation systems describing the energization of a transformer, using appropriate numerical method for the correct simulations of equations. In general, for simulating electrical systems, the priority is that the numerical method is A-stable, has a high accuracy, and successfully solves a (very) stiff system.

In regard to these requirements, the conclusion is that methods BDF1 and -2 and NDF1 and -2 are A- and L-stable. They are generally well suited for simulations of nonlinear electrical systems. This paper realizes an algorithm for the solution of (7) based on the usages of the NDF2 or BDF2 methods. These methods were selected because of the fact that they have the same truncation errors of order At3, like a trapezoidal method, and they are A- and L-stable. On the other hand, it is shown that it is possible to combine these methods using the widely used trapezoidal method for simulating electrical systems.

4. Inrush Current

Measurements and Simulations

Three developed single-phase transformer models were tested on an example of a transformer inrush current transient.

The system and transformer parameters are

(i) system voltage: e(t) = 328 cos(o>£ - 39°);

(ii) nominal transformer power: Str = 300 VA;

(iii) primary/secondary voltage: Up = 220/24 V;

(iv) short circuit voltage: uk% = 5.45%;

(v) primary winding resistance: Rp = 3.75 O;

(vi) primary winding leakage inductance: Lp = 2.54 mH;

(vii) core loss resistance: Rm = 3826 O.

The nonlinear magnetizing curve with a major hysteresis loop is shown in Figure 9. A nonlinear curve of an iron core losses resistor is shown in Figure 10.

Figures 11 and 12 present the results of measurement of transformer inrush and steady-state current. In addition, harmonic analysis of the transformer inrush current for the first four harmonic components is realized, Figure 13.

Figure 14 shows the simulation results of the developed algorithm for Model 3 realized using the classical trapezoidal method with integration step At = 100 ^sec.

Figure 14 clearly shows the existence of numerical oscillations during the trapezoidal method's usage. Otherwise, any simulations of the developed models (Model 1, Model 2, or Model 3), using trapezoidal method, clearly showed the existence of unwanted numerical oscillations.

In order to investigate the causes of numerical oscillations during simulation of maximum and minimum eigenvalues, stiffness ratios and stiffness indexes were calculated at each integration step. Figure 15 shows the results of the calculated maximum and minimum eigenvalues during the simulation times.

x10 -1.5

Min eigenvalue

Time (s) Model 1

Max eigenvalue

T-i—TT—T7—

x10 -0.5

Min eigenvalue

-4.5 0.14 0

Time (s) Model 2

Max eigenvalue

-100 0.14 0

x10 -4.5

Min eigenvalue

-4.9 0.14 0

-100 0.14 0

Time (s) Time (s) Time (s)

Model 1 Model 2 - Model 3

Figure 15: Eigenvalue propagation during the simulation time: all three models.

Time (s) Model 3

Max eigenvalue

i"\i"\f\r\!"\!

Inrush current (measurement and simulation)

Steady-state current (measurement and simulation)

Measured Model 1

0.06 0.08 Time (s)

Model 2 - Model 3

Figure 16: Measured and simulated transformer inrush current, NDF2 method, At = 100 ^sec.

Measured Model 1

0.12 0.13 0.14 Time (s)

Model 2 - Model 3

Figure 17: Measured and simulated transformer steady-state current, NDF2 method, At = 100 ^sec.

Table 1 shows the borders of the analyzed stiffness parameters. It is obvious that they are very stiff systems for all three cases.

It is interesting that the extreme values of parameters rose in the unsaturated regions of the magnetizing (hysteresis) curve, and the extreme values of parameter ^ k

rose in the saturated region of this curve. Also, it is interesting that the time function of parameter ^ k is an analog to the inrush current waveform in all three cases.

In view of the mentioned reasons, the BDF2 and NDF2 methods were used for simulation of transformer inrush current. Figures 16 and 17 show the results when comparing

34 Peak number

Model 1 -<— Model 2 -e- Model 3

Figure 18: Relative error of peak values of inrush current.

3.5 3 2.5 2 1.5 1

3.5 3 2.5 2 1.5 1

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Time (s)

0.02 0.04 0.06 0.08 0.1 0.12 0.14 Time (s)

3.5 3 2.5 2 1.5 1

3.5 3 2.5 2 1.5 1

0.02 0.04 0.06 0.08 0.1 0.12 0.14 Time (s)

0.02 0.04 0.06 0.08 0.1 0.12 0.14 Time (s)

Measured Model 1

Model 2 Model 3

Measured Model 1

Model 2 Model 3

Figure 19: Harmonic content of measured and simulated inrush currents.

Table 1: Stiffness ratios and stiffness indexes of the system.

Model 1 1.508 • 1.965 • ' 106 < Cj.fc < 1.596 • 106 104 < ^ < 4.949 • 104

Model 2 0.964 • 1.256 • ■ 106 < (yt < 3.775 • 106 104 < < 1.183 • 107

Model 3 4.521 • 4.622 • ' 106 < Cj.fc < 4.784 • 106 104 < < 1.566 • 107

Table 2: Peak error of steady-state current.

Model Model1 Model 2 Model 3

Peak error [%] 13.18% 11.98% 4.56%

the measured and simulated inrush current and steady-state transformer current by usage of the NDF2 method with an integration step of Af = 100 ^sec.

Figure 18 shows the peak values for the inrush currents errors for all three models. It is necessary to note that the maximum error when using Models 1 and 2 was 5.29%, whilst when using Model 3, the maximum error was 4.30%.

Table 2 shows the error in peak value for the transformer steady-state current. Minimum error was in Model 3.

In addition, a comparison of the harmonic contents between the measured and simulated inrush current is realized, Figure 19. This figure shows that the best results were for Model 3, whilst Models 1 and 2 provided mostly identical simulation results.

On the basis of all the aforementioned, it is possible to suggest the hybrid numerical method as a linear combination of the traditional trapezoidal method and the proposed BDF2 (NDF2) method. In [38] a hybrid method was proposed by Tadeusiewicz and Halgas in 2005 that presented a linear combination of the trapezoidal and BDF2 method in the form

( 0 \ 0 Ai

X„+1 = (1 + 3jx„ + (1 - 0) —F (X„, o

+ (l + 0) A~F(X„+I,i„+i).

In expression (14), 0 = 0 leads to the trapezoidal method, whereas 0=1 leads to the BDF2 method.

In the same way, we have proposed a hybrid method that represents a linear combination of the trapezoidal and NDF2 methods in the form

/ 0 \ 30 Ai

X„+1 = (1 + X„ - -X„-1 + (1 - 0) y F (X„, O

+ (1 + 5)|F (*»+i> f„+i) + 10xio+]i-

In the expression (15), 0 = 0 leads to the trapezoidal method, whereas 0 = 1 leads to the NDF2 method. This method overcomes all the eventual problems during the usage of the standard trapezoidal method.

5. Conclusions

This paper investigates the different models and numerical methods that could be implemented for the simulations of inrush current in a single-phase transformer. Various core models of a transformer are developed in the paper: a model with a nonlinear inductor and a linear resistor, a model with a nonlinear inductor and a nonlinear resistor, and a model with a hysteresis inductor and a linear resistor. In addition, the paper investigates eigenvalues propagations during simulation time for all three models. The defined stiffness parameters, the stiffness ratios and the stiffness indexes, indicated that very stiff systems should be solved for all three developed models. The state-space equation system is solved by using proposed A- and L-stable numerical differentiation formulas (NDF2) (BDF2) method. The proposed method has the same order as a trapezoidal method; however this method overcomes the main drawback of the trapezoidal method (numerical oscillations). In the case where the numerical oscillations are generated by use of the trapezoidal method, the NDF2 (BDF2) method efficiently suppresses them. Comparisons between the measured and simulated inrush and steady-state currents of a transformer are conducted for all three proposed models. The best simulation results were obtained by using the developed Model 3, which incorporates hysteresis inductor. The developed Models 1 and 2 gave almost identical simulation results. The realized nonlinear inductor, nonlinear resistor, and hysteresis model can be incorporated within the EMTP-type of programs by using a combination of the existing trapezoidal and proposed NDF (BDF) method of order 2. In this way, the suggested numerical method could be used in the simulation of electrical networks with nonlinear elements. The future work will investigate the propagation of eigenvalues of a system in simulations of inrush currents in a three-phase power transformer.

References

[1] CIGRE Working Group SC 33.02, "Guidelines for representation of network elements when calculating transients," CIGRE Brochure, vol. 39, pp. 1-29,1990.

[2] D. A. N. Jacobson, P. W. Lehn, and R. W. Menzies, "Stability domain calculations of period-1 ferroresonance in a nonlinear resonant circuit," IEEE Transactions on Power Delivery, vol. 17, no. 3, pp. 865-871, 2002.

[3] G. Stumberger, S. Seme, B. Stumberger, B. Polajzer, and D. Dolinar, "Determining magnetically nonlinear characteristics of transformers and iron core inductors by differential evolution," IEEE Transactions on Magnetics, vol. 44, no. 6, pp. 15701573, 2008.

[4] C. Perez-Rojas, "Fitting saturation and hysteresis via arctangent functions," IEEE Power Engineering Review, vol. 20, no. 11, pp. 55-57, 2000.

[5] X. Zhang, "An efficient algorithm for transient calculation of the electric networks including magnetizing branches," Mathematical Problems in Engineering, vol. 2011, Article ID 479626, 11 pages, 2011.

[6] A. Ben-Tal, V. Kirk, and G. Wake, "Banded chaos in power systems," IEEE Transactions on Power Delivery, vol. 16, no. 1, pp. 105-110, 2001.

[7] H. W. Dommel, Electromagnetic Transients Program, Reference Manual, (EMTP Theory Book), Bonneville Power Administration, Portland, Ore, USA, 1986.

[8] C. E. Lin, C. L. Cheng, C. L. Huang, and J. C. Yeh, "Investigation of magnetizing inrush current in transformers—part I: numerical simulation," IEEE Transactions on Power Delivery, vol. 8, no.

1, pp. 246-254,1993.

[9] E. H. Badawy and R. D. Youssef, "Representation oftransformer saturation," Electric Power Systems Research, vol. 6, no. 4, pp. 301-304,1983.

[10] S. Seker, T. C. Akinci, and S. Taskin, "Spectral and statistical analysis for ferroresonance phenomenon in electric power systems," Electrical Engineering, vol. 94, no. 2, pp. 117-124, 2012.

[11] I. Clavería, M. García-Gracia, M. A. García, and L. Montanes, "A time domain small transformer model under sinusoidal and non-sinusoidal supply voltage," European Transactions on Electrical Power, vol. 15, no. 4, pp. 311-323, 2005.

[12] J. R. Lucas, P. G. McLaren, W. W. L. Keerthipala, and R. P. Jayasinghe, "Improved simulation models for current and voltage transformers in relay studies," IEEE Transactions on Power Delivery, vol. 7, no. 1, pp. 152-159,1992.

[13] Y. Baghzouz and X. D. Gong, "Voltage-dependent model for teaching transformer core nonlinearity," IEEE Transactions on Power Systems, vol. 8, no. 2, pp. 746-752,1993.

[14] T. Tran-Quoc and L. Pierrat, "Efficient non linear transformer model and its application to ferroresonance study," IEEE Transactions on Magnetics, vol. 31, no. 3, pp. 2060-2063,1995.

[15] A. D. Theocharis, J. Milias-Argitis, and T. Zacharias, "Singlephase transformer model including magnetic hysteresis and eddy currents," Electrical Engineering, vol. 90, no. 3, pp. 229-241, 2008.

[16] J. Faiz and S. Saffari, "Inrush current modeling in a single-phase transformer," IEEE Transactions on Magnetics, vol. 46, no. 2, pp. 578-581, 2010.

[17] H. Akçay and D. G. Ece, "Modeling of hysteresis and power losses in transformer laminations," IEEE Transactions on Power Delivery, vol. 18, no. 2, pp. 487-492, 2003.

[18] S.-R. Huang, S. C. Chung, B.-N. Chen, and Y.-H. Chen, "A harmonic model for the nonlinearities of single-phase transformer with describing functions," IEEE Transactions on Power Delivery, vol. 18, no. 3, pp. 815-820, 2003.

[19] W. L. A. Neves and H. W. Dommel, "On modelling iron core nonlinearities," IEEE Transactions on Power Systems, vol. 8, no.

2, pp. 417-422, 1993.

[20] A. Tokic, I. Uglesic, V. Milardic, and G. Stumberger, "Conversion of RMS to instantaneous saturation curve: inrush current and ferroresonance cases," in Proceedings of the 14th International IGTE Symposium on Numerical Field Calculation in Electrical Engineering, pp. 1-4, Graz, Austria, 2010.

[21] I. D. Mayergoyz, Mathematical Models of Hysteresis and Their Applications, Academic Press, Elsevier, New York, NY, USA, 2003.

[22] F. C. Trutt, E. A. Erdelyi, and R. E. Hopkins, "Representation of the magnetization characteristic of DC machines for computer use," IEEE Transactions on Power Apparatus and Systems, vol. 87, no. 3, pp. 665-669,1968.

[23] O. Karlqvist, Calculation of the Magnetic Field in the Ferromagnetic Layer of a Magnetic Drum, vol. 86 of Transactions of the Royal Institute of Technology, Elanders boktr., Stockholm, Sweden, 1954.

[24] J. Takacs, "A phenomenological mathematical model of hysteresis," COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, vol. 20, no. 4, pp. 1002-1014, 2001.

[25] J. Takacs, "The Everett integral and its analytical approximation," in Advanced Magnetic Materials, pp. 203-230, Intech Publication, 2012.

[26] A. Tokic, V. Madzarevic, and I. Uglesic, "Numerical calculations of three-phase transformer transients," IEEE Transactions on Power Delivery, vol. 20, no. 4, pp. 2493-2500, 2005.

[27] A. Tokic and I. Uglesic, "Elimination of overshooting effects and suppression of numerical oscillations in transformer transient calculations," IEEE Transactions on Power Delivery, vol. 23, no. 1, pp. 243-251, 2008.

[28] A. Tokic, V. Madzarevic, and I. Uglesic, "Hysteresis model in transient simulation algorithm based on BDF numerical method," in Proceedings of the IEEE Power Tech Conference, pp. 1-6, St. Petersburg, Russia, 2005.

[29] E. Hairer and G. Wanner, Solving Ordinary Differential Equations. II: Stiff and Differential-Algebraic Problems, vol. 14 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1991.

[30] C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, NJ, USA, 1971.

[31] R. W. Klopfenstein, "Numerical differentiation formulas for stiff systems of ordinary differential equations," RCA Review, vol. 32, no. 3, pp. 447-462,1971.

[32] L. F. Shampine and M. W. Reichelt, "The MATLAB ODE suite," SIAM Journal on Scientific Computing, vol. 18, no. 1, pp. 1-22, 1997.

[33] J. R. Cash, "On the integration of stiff systems of O.D.E.s using extended backward differentiation formulae," Numerische Mathematik, vol. 34, no. 3, pp. 235-246,1980.

[34] C. Fredebeul, "A-BDF: a generalization of the backward differentiation formulae," SIAM Journal on Numerical Analysis, vol. 35, no. 5, pp. 1917-1938,1998.

[35] J. E. Frank and P. J. van der Houwen, "Diagonalizable extended backward differentiation formulas," BIT Numerical Mathematics, vol. 40, no. 3, pp. 497-512, 2000.

[36] E. Alberdi and J. J. Anza, "A predictor modification to the EBDF method for Stiff systems," Journal of Computational Mathematics, vol. 29, no. 2, pp. 199-214, 2011.

[37] S. A. M. Yatim, Z. B. Ibrahim, K. I. Othman, and M. B. Suleiman, "A quantitative comparison of numerical method for solving stiff ordinary differential equations," Mathematical Problems in Engineering, vol. 2011, Article ID 193691, 12 pages, 2011.

[38] M. Tadeusiewicz and S. Halgas, "Transient analysis of nonlinear dynamic circuits using a numerical-integration method," COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, vol. 24, no. 2, pp. 707719, 2005.

Copyright of Mathematical Problems in Engineering is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.