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

Authors: Amir Tokić; Ivo Uglešić; Gorazd Štumberger
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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 usingA- 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.


---

## Body

## 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 Π-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 Π equivalent model.Figure 1
StandardT equivalent transformer circuit model.Different types ofT 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 resistorRm and nonlinear inductor Lm,(b)
nonlinear resistorRm and nonlinear inductor Lm,(c)
linear resistorRm and nonlinear hysteresis inductor Lm.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, [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 Takâcs in [24] and refined by Takâcs 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.


### 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 Figure2.Figure 2
Nonlinear curve of core inductor.When this curve is piecewise linear, the results are the input vectors magnetizing currentsILm=[is1,is2,…,isN]T and magnetic fluxes ΦLm=[ϕs1,ϕs2,…,ϕsN]T, whereas N is the total number of magnetizing curve regions. The magnetizing current for the kth linear region of the magnetizing curve, k=1,2,…,(N-1), is calculated using the equation
(1)imk=1Lmkϕ+sign(ϕ)Isk=1Lmkϕ+Sϕk,
where Lmk=(ϕsk+1-ϕsk)/(isk+1-isk)andIsk=isk-(1/Lmk)ϕsk.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.(a) Nonlinear inductor and (b) equivalent model.
(a)
(b)


### 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 Figure4.Figure 4
Nonlinear curve of core resistor.In this case, the input vectors are the resistor currentsIRm=[ir1,ir2,…,irM]T and resistor voltages URm=[ur1,ur2,…,urM]T, where M is the total number of resistor curve regions. The resistor current for the kth linear region of the resistor curve, k=1,2,…,(M-1), is calculated using the equation
(2)iRk=1RmkuR+sign(uR)Irk=1RmkuR+Srk,
where Rmk=(urk+1-urk)/(irk+1-irk)  and  Irk=irk-(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, Figure5.(a) Nonlinear resistor and (b) equivalent model.
(a)
(b)


### 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 (Figure6) (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.Figure 6
Major and minor hysteresis loops.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 currentsILh=[ih1,ih2,…,ihP]T  and magnetic fluxes  ΦLh=[ϕh1,ϕh2,…,ϕhP]T, where P is the total number of the major hysteresis loops.The expression for the hysteresis current of thekth piecewise region in terms of the main flux  ϕm, k=1,2,…,(P-1), Figure 6, can be stipulated as
(3)ihk=1Lmkϕm+sign(Δϕm)Ihk,
where Δϕm=ϕm(tj)-ϕm(tj-1), j=1,2,3,…,Jmax, Lmk=(ϕhk+1-ϕhk)/(ihk+1-ihk),  and  Ihk=ihk-(1/Lmk)ϕhk.In the general case, if it is assumed that the distance between the operating point and the major loop is a linear function of the magnetic flux, for both the ascending and descending trajectory, the next equation is valid:(4)d=μkϕ+ηk=sign(Δϕ)(ϕ-ϕm).The unknown coefficientsμk and ηk 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
(5)ihk=1Lmk-ΔLmkϕ+sign(Δϕ)·(Ihk-ΔIk),
noting that ΔLmk=(sign(Δϕ)μk/(sign(Δϕ)μk-1))Lmk and ΔIk=(1/Lmk)ηk.In this way the hysteresis inductorLh in Figure 3(a) can be modeled by introducing a linear inductor Lhk in parallel connection with a current source  Shk, both functionally dependent on the actual linear region k, Figure 7 as follows:
(6)Lhk=Lmk-ΔLmk,Shk=sign(Δϕ)·(Ihk-ΔIk).(a) Hysteresis inductor and (b) equivalent model.
(a)
(b)It is noted that all three developed models has a special routine for eliminating the possible overshooting effect [26–28].


## 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 Figure2.Figure 2
Nonlinear curve of core inductor.When this curve is piecewise linear, the results are the input vectors magnetizing currentsILm=[is1,is2,…,isN]T and magnetic fluxes ΦLm=[ϕs1,ϕs2,…,ϕsN]T, whereas N is the total number of magnetizing curve regions. The magnetizing current for the kth linear region of the magnetizing curve, k=1,2,…,(N-1), is calculated using the equation
(1)imk=1Lmkϕ+sign(ϕ)Isk=1Lmkϕ+Sϕk,
where Lmk=(ϕsk+1-ϕsk)/(isk+1-isk)andIsk=isk-(1/Lmk)ϕsk.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.(a) Nonlinear inductor and (b) equivalent model.
(a)
(b)


## 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 Figure4.Figure 4
Nonlinear curve of core resistor.In this case, the input vectors are the resistor currentsIRm=[ir1,ir2,…,irM]T and resistor voltages URm=[ur1,ur2,…,urM]T, where M is the total number of resistor curve regions. The resistor current for the kth linear region of the resistor curve, k=1,2,…,(M-1), is calculated using the equation
(2)iRk=1RmkuR+sign(uR)Irk=1RmkuR+Srk,
where Rmk=(urk+1-urk)/(irk+1-irk)  and  Irk=irk-(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, Figure5.(a) Nonlinear resistor and (b) equivalent model.
(a)
(b)


## 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 (Figure6) (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.Figure 6
Major and minor hysteresis loops.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 currentsILh=[ih1,ih2,…,ihP]T  and magnetic fluxes  ΦLh=[ϕh1,ϕh2,…,ϕhP]T, where P is the total number of the major hysteresis loops.The expression for the hysteresis current of thekth piecewise region in terms of the main flux  ϕm, k=1,2,…,(P-1), Figure 6, can be stipulated as
(3)ihk=1Lmkϕm+sign(Δϕm)Ihk,
where Δϕm=ϕm(tj)-ϕm(tj-1), j=1,2,3,…,Jmax, Lmk=(ϕhk+1-ϕhk)/(ihk+1-ihk),  and  Ihk=ihk-(1/Lmk)ϕhk.In the general case, if it is assumed that the distance between the operating point and the major loop is a linear function of the magnetic flux, for both the ascending and descending trajectory, the next equation is valid:(4)d=μkϕ+ηk=sign(Δϕ)(ϕ-ϕm).The unknown coefficientsμk and ηk 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
(5)ihk=1Lmk-ΔLmkϕ+sign(Δϕ)·(Ihk-ΔIk),
noting that ΔLmk=(sign(Δϕ)μk/(sign(Δϕ)μk-1))Lmk and ΔIk=(1/Lmk)ηk.In this way the hysteresis inductorLh in Figure 3(a) can be modeled by introducing a linear inductor Lhk in parallel connection with a current source  Shk, both functionally dependent on the actual linear region k, Figure 7 as follows:
(6)Lhk=Lmk-ΔLmk,Shk=sign(Δϕ)·(Ihk-ΔIk).(a) Hysteresis inductor and (b) equivalent model.
(a)
(b)It is noted that all three developed models has a special routine for eliminating the possible overshooting effect [26–28].


## 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. 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, Figure8.Figure 8
Equivalent transformer model for inrush current simulations.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:
(7)dX(t)dt=Aj,kX(t)+Bj,kU(t)=F(X,t),
where the state-space vector and input vector are:
(8)X(t)=[iLpϕ]T,U(t)=[e(t)SϕkSrj]T.System matrixes are, respectively as follows:(9)Aj,k=[-Rp+RmjLpRmjLpLmkRmj-RmjLmk],Bj,k=[1Lp0RmjLp-RmjRmjLp-Rmj]T.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ξj,k and stiffness index ζj,k for all state matrixes Aj,k.:
(10)ξj,k=maxi|Re(λj,k(i))|mini|Re(λj,k(i))|,ζj,k=maxi|Re(λj,k(i))|,
where λj,k(i), i=1,2,…,dim(Aj,k) are eigenvalues of all state matrixes Aj,k.. In the case for ξj,k≫1 it is a stiff system, and in the case for ζj,k→∞ it is a very stiff system.The state equations (7) are traditionally solved using the implicit A-stable, second-order trapezoidal method [7]:
(11)Xn+1=Xn+Δt2[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 Aj,k 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 (BDFp) methods of pth order, p≤5, as introduced by Gear in 1971, have the following form [29, 30]:
(12)∑m=1p1m∇mXn+1=ΔtF(Xn+1,tn+1).The order of the truncation errors for the BDFp methods is Δtp+1 [29, 30]. The BDFp methods are A(α)- and L-stable, with stability angles ranging between 90°, 90°, 86°, 73°, and 51° [25]. The L-stability properties of the BDFp methods lead to the suppression of numerical oscillations, 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]:
(13)∑m=1p1m∇mXn+1=ΔtF(Xn+1,tn+1)FFfffffffffffffffFFF+κγp(Xn+1-Xn+1[0]),
where the coefficients are γp=∑m=1p(1/m), the predicted value is Xn+1[0]=∑m=0p∇mXn, and parameter κ 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(α)- 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 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 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 isA-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 areA- 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 Δt3, 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)=328cos(ωt-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 Ω;(vi)
primary winding leakage inductance:Lp = 2.54 mH;(vii)
core loss resistance:Rm = 3826 Ω.The nonlinear magnetizing curve with a major hysteresis loop is shown in Figure9. A nonlinear curve of an iron core losses resistor is shown in Figure 10.Figure 9
Magnetizing curve and hysteresis loop.Figure 10
Nonlinear core loss resistor curve.Figures11 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 11
Measured transformer inrush current.Figure 12
Measured transformer steady-state current.Figure 13
Harmonic content of measured inrush current.Figure14 shows the simulation results of the developed algorithm for Model 3 realized using the classical trapezoidal method with integration step Δt=100 μsec.Figure 14
Inrush current simulation, Model 3, trapezoidal method,Δt=100 μsec.Figure14 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. Figure15 shows the results of the calculated maximum and minimum eigenvalues during the simulation times.Figure 15
Eigenvalue propagation during the simulation time: all three models.Table1 shows the borders of the analyzed stiffness parameters. It is obvious that they are very stiff systems for all three cases.Table 1
Stiffness ratios and stiffness indexes of the system.
Model 1
1.508
·
1
0
6
≤
ζ
j
,
k
≤
1.596
·
1
0
6
1.965
·
1
0
4
≤
ξ
j
,
k
≤
4.949
·
1
0
4
Model 2
0.964
·
1
0
6
≤
ζ
j
,
k
≤
3.775
·
1
0
6
1.256
·
1
0
4
≤
ξ
j
,
k
≤
1.183
·
1
0
7
Model 3
4.521
·
1
0
6
≤
ζ
j
,
k
≤
4.784
·
1
0
6
4.622
·
1
0
4
≤
ξ
j
,
k
≤
1.566
·
1
0
7It is interesting that the extreme values of parametersξj,k rose in the unsaturated regions of the magnetizing (hysteresis) curve, and the extreme values of parameter ζj,k rose in the saturated region of this curve. Also, it is interesting that the time function of parameter ζj,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. Figures16 and 17 show the results when comparing the measured and simulated inrush current and steady-state transformer current by usage of the NDF2 method with an integration step of Δt=100 μsec.Figure 16
Measured and simulated transformer inrush current, NDF2 method,Δt=100 μsec.Figure 17
Measured and simulated transformer steady-state current, NDF2 method,Δt=100 μsec.Figure18 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%.Figure 18
Relative error of peak values of inrush current.Table2 shows the error in peak value for the transformer steady-state current. Minimum error was in Model 3.Table 2
Peak error of steady-state current.
Model
Model 1
Model 2
Model 3
Peak error[%]
13.18%
11.98%
4.56%In addition, a comparison of the harmonic contents between the measured and simulated inrush current is realized, Figure19. This figure shows that the best results were for Model 3, whilst Models 1 and 2 provided mostly identical simulation results.Figure 19
Harmonic content of measured and simulated inrush currents.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
(14)Xn+1=(1+θ3)Xn-θ3Xn-1+(1-θ)Δt2F(Xn,tn)+(1+θ3)Δt2F(Xn+1,tn+1).In expression (14), θ=0 leads to the trapezoidal method, whereas θ=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(15)Xn+1=(1+θ5)Xn-3θ10Xn-1+(1-θ)Δt2F(Xn,tn)+(1+θ5)Δt2F(Xn+1,tn+1)+θ10Xn+1[0].In the expression (15), θ=0 leads to the trapezoidal method, whereas θ=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 proposedA- 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.



---

*Source: 215647-2013-09-11.xml*