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 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 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.
Different types of 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, and represent the serial winding resistances that generally include the Joule and eddy current losses in the windings, whilst and represent the serial leakage inductances divided amongst primary and secondary windings. The shunt branch parameters describe iron core behavior, where the resistance represents the core eddy current losses, whilst inductance 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 and nonlinear inductor ,(b)nonlinear resistor and nonlinear inductor ,(c)linear resistor and nonlinear hysteresis inductor .
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 Figure 2.
When this curve is piecewise linear, the results are the input vectors magnetizing currents and magnetic fluxes , whereas is the total number of magnetizing curve regions. The magnetizing current for the th linear region of the magnetizing curve, , is calculated using the equation where and.
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 . At the same time, contrary to EMTP-based equivalent models, the inductor model is not dependent on the integration step.
(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 Figure 4.
In this case, the input vectors are the resistor currents and resistor voltages , where is the total number of resistor curve regions. The resistor current for the th linear region of the resistor curve, , is calculated using the equation where and .
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.
(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 (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 and magnetic fluxes , where is the total number of the major hysteresis loops.
The expression for the hysteresis current of the th piecewise region in terms of the main flux , , Figure 6, can be stipulated as where , , , and .
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:
The unknown coefficients and are determined on condition that the two successive reversal points and lie on this line.
Based on relations (3)-(4) it can be proved that the expression for the hysteresis current of the th linear region that describes the trajectory of the actual operating point is noting that and .
In this way the hysteresis inductor in Figure 3(a) can be modeled by introducing a linear inductor in parallel connection with a current source , both functionally dependent on the actual linear region , Figure 7 as follows:
(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, 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: where the state-space vector and input vector are:
System matrixes are, respectively as follows:
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 and stiffness index for all state matrixes .: where , are eigenvalues of all state matrixes .. In the case for it is a stiff system, and in the case for it is a very stiff system.
The state equations (7) are traditionally solved using the implicit -stable, second-order trapezoidal method [7]:
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 -stable, such that, for state matrixes 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) methods of th order, , as introduced by Gear in 1971, have the following form [29, 30]:
The order of the truncation errors for the BDF methods is [29, 30]. The BDF methods are - and -stable, with stability angles ranging between 90°, 90°, 86°, 73°, and 51° [25]. The -stability properties of the BDF 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 (NDF), as suggested by Klopfenstein in 1971, are a modification of the BDF’s methods with the next relations [31, 32]: where the coefficients are , the predicted value is , 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 - and -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 -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 is -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 - and -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 , like a trapezoidal method, and they are - and -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: ;(ii)nominal transformer power: = 300 VA;(iii)primary/secondary voltage: = 220/24 V;(iv)short circuit voltage: % = 5.45%;(v)primary winding resistance: = 3.75 Ω;(vi)primary winding leakage inductance: = 2.54 mH;(vii)core loss resistance: = 3826 Ω.
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 μ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.
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 rose in the saturated region of this curve. Also, it is interesting that the time function of parameter 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 the measured and simulated inrush current and steady-state transformer current by usage of the NDF2 method with an integration step of μ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
In expression (14), leads to the trapezoidal method, whereas 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
In the expression (15), leads to the trapezoidal method, whereas 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 - and -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.