An Efficient Algorithm for Transient Calculation of the Electric Networks Including Magnetizing Branches

A simplified model of magnetic saturation characteristics is proposed in this paper for transient calculation of the electric networks including magnetizing branches. The model represents the magnetic saturation characteristics by the continuous function instead of the piecewise linear approximation. Based on the proposed model, an efficient transient algorithm is developed. The nonlinear differential equations describing the transient behavior of the magnetizing branches are solved by the semiexplicit Runge-Kutta method, in which noniterative computations are involved. The transient calculation for the remaining linear network is performed in terms of the solution to the magnetizing branches. A comparison is made between calculated and experimental results to confirm the validity of the algorithm.


Introduction
Electric power systems contain a series of ferromagnetic components, such as power transformers, potential transformers, and shunt reactors.Owing to the presence of the ironcores, the ferromagnetic components have nonlinear magnetizing characteristics.When a switching or fault operation occurs in the vicinity of the ferromagnetic components, their magnetic saturation characteristics may exert a significant influence on the transient phenomena.In such a situation, it is necessary to take account of the magnetic saturation characteristics in transient calculation.The magnetic saturation characteristics are usually described by the flux linkagemagnetizing current curves, namely, the saturation curves.The piecewise linear approach has been widely employed for simulating the saturation curves in the transient solution of magnetizing branches 1-5 .In fact, this approach can only give a roughly approximate treatment to the nonlinearity of the saturation curves.When a saturation curve is approximated in piecewise linear form, some illegal overshoot conditions are frequently caused at the turning points of piecewise linear inductances, which may give rise to a numerical error for the transient calculation 4 .The compensation method has also been used to solve the transient phenomena associating with nonlinear elements 4, 5 .However, it leads to an increase in computational burden due to its iterative solution.In view of these drawbacks, this paper proposes a simplified model of magnetic saturation characteristics for transient solution of magnetizing branches.The saturation curves are represented by the continuous function instead of the piecewise linearization.On the basis of the model, an efficient algorithm is developed for transient calculation.In the algorithm, the electric network is divided into two parts, that is, the magnetizing branches and the remaining linear network.The nonlinear differential equations are formed for describing the transient behavior of the magnetizing branches.The semiexplicit Runge-Kutta method is used to numerically solve these equations, which can avoid the iterative computations.The transient calculation for the remaining linear network can be performed subsequent to the solution to the magnetizing branches.In order to confirm the validity of the algorithm, the calculated results are compared with the experimental ones.

Modeling of Magnetizing Branches
From the point of view of the practical application in transient calculation, a magnetizing branch can be simplified as a nonlinear inductance in parallel with an appropriate core loss resistance R m 6-9 , as shown in Figure 1.The constant resistance R m can be determined from the manufacturer's data on the core loss of ferromagnetic components versus applied voltage 8, 9 .The nonlinear inductance is characterized by the saturation curve, as shown in Figure 2.
In this study, the saturation curve is formulated as where λ is flux linkage, i is magnetizing current, and A, B, W are the constants which can be determined from the corresponding measured data by using the least squares curve fitting method.The dynamic inductance necessary to transient calculation can be derived from 2.1 In the unsaturation region, the small magnetizing current yields a large linear inductance In the saturation region, the magnetizing current becomes large and the dynamic inductance is reduced to It is clear from 2.3 and 2.4 that 2.1 can give a reasonable description to the nonlinear saturation behavior of practical ferromagnetic components, such as power transformers and potential transformers.

Algorithm for Transient Calculation
An electric network including N magnetizing branches is depicted in Figure 3, where the N core loss resistances R m have been incorporated into the linear network.Let λ and i denote the flux linkage and magnetizing current vectors, respectively, that is, λ  a Thevenin equivalent procedure is made for Figure 3 1, 4 .By removing the N nonlinear inductances from Figure 3, the opened-circuit voltage vector u 0 u 0 u 10 , u 20 , . . ., u N0 T and input impedance matrix Z 0 N × N associated with the N ports can be conveniently calculated by the method given in 4, 5 , as shown in Figure 4.After obtaining u 0 and Z 0 , Figure 3 is simplified as a Thevenin equivalent network expressed in matrix form, as shown in Figure 5.The matrix differential equation can be given for the N nonlinear inductances dλ dt u 0 − Z 0 i.

3.1
The above derivative can be changed into The nonlinear dynamic inductance matrix L d i is given as where diag denotes the diagonal matrix and the N nonlinear inductances are considered to be uncoupling for the real reason.Each diagonal element of the matrix L d i can be determined by 2.2 .Substituting 3.2 into 3.1 leads to Equation 3.5 represents N nonlinear differential equations and may be a system of stiff differential equations when modeling the real data of the electric networks.This kind of problem is apparent in numerical solution of power transformer transients.The stiffness of differential equations makes it difficult for the explicit difference scheme to solve these equations successfully.The explicit difference scheme, applied to stiff differential equations, is numerically unstable because this difference scheme results in an increase of truncation error 10 .The implicit difference scheme, although numerically stable, does increase the solution time for the magnetizing branches due to its incorporating of the iterative computations into the numerical solution process.Obviously, a suitable solving technique for 3.5 should be numerically stable and avoid the iterative computations.The semiexplicit Runge-Kutta method can accomplish this purpose.A generalized semiexplicit Runge-Kutta difference scheme is given as follows 11 :

3.6
where A is the Jacobian matrix, A i ∂F i /∂i.Through a quite tedious calculation, it is possible to expand i n 1 − i n in 3.6 as a power series in Δt and compare this with the Taylor's series.The detail was given in 12 .In order to ensure the correspondence between the early terms of both the series, the following relationships should be satisfied for a truncation error 0 Δt 4 :

3.7
These relationships involve six adjustable constants, k 1 , k 2 , a 1 , a 2 , b 1 , and c 1 .Choosing the first two constants k 1 0.75 and k 2 0.25 13 , the remaining four constants can be determined as

3.8
Thus, 3.6 is simplified as a two-stage and third-order difference scheme where 1 is unit matrix.Once the magnetizing current vector i is obtained by using the difference formula 3.9 to solve 3.5 , the transient calculation of the linear network in Figure 3 i r1 λ r1 λ t3

Transmission lines Coupled branches
Magnetizing branches at receiving end Magnetizing branches at sending end  can be performed easily.By replacing the N nonlinear inductances with the N current sources equal to the N known elements of the magnetizing current vector i, as shown in Figure 6, Figure 3 is changed into a completely linear network.The transient responses can be calculated for the linear network by using the algorithm presented in 1, 4, 5 .

Case Study
An electric network with two groups of magnetizing branches at both ends of transmission lines is shown in Figure 7, where the group at the sending end corresponds to a transformer and the other at the receiving end corresponds to a reactor.The core loss resistance of the reactor is ignored 14 .The network data are given in the following:  For the sake of convenience, Figure 7 is depicted as a single phase network taking matrix form, as shown in Figure 8.The R-L coupled branches and the transmission lines are further represented as their respective transient equivalent circuits, as shown in Figure 9.
In the transient equivalent circuits, R RL and I RL t are the equivalent resistance matrix and current source vector of the R-L coupled branches, and Z, I ls t and I lr t are the phase surge impedance matrix and historical current source vectors of the transmission lines, respectively.The special formulas for calculating these matrices and vectors have been derived in 5, 15 .By combining Figure 9 with Figure 8, an equivalent calculation network is given in Figure 10.
With the terminals of the transformer and reactor magnetizing branches open-circuited in Figure 10, two Thevenin equivalent subnetworks can be found, and then the transient calculation can be performed by following the procedure stated above.The calculated voltage waveforms at the receiving end of transmission lines are shown in Figure 11, where the corresponding experimental waveforms obtained from the laboratory-scale model network are simultaneously given for comparison.It can be seen from Figure 11 that a better consistency appears between calculated and experimental results.In addition, the receiving end voltages are calculated by the traditional piecewise linear approach.The nonlinear saturation curve is simulated as piecewise linear inductances with 3 slopes.The calculated result shows that sudden jumps occur on the calculated voltage waveforms and make them unsmoothed, as illustrated in Figure 12.However, such a numerical error is not found on the voltage waveforms calculated by the proposed algorithm.The voltage peak values calculated by the proposed algorithm are closer to the experimental values than these by the piecewise linear approach.

Conclusions
An approach to modeling the magnetic saturation characteristics of magnetizing branches has been presented.It uses a continuous function to represent the saturation curves instead of the piecewise linearization and is capable of characterizing the saturation nonlinearity of practical ferromagnetic components.By implementing this approach into transient calculation, an efficient algorithm has been developed.The transient behavior of the magnetizing branches is described by a set of nonlinear differential equations, and the semiexplicit Runge-Kutta method is employed to solve these equations so that the iterative computations can be avoided.Then, the complete network solution can be obtained in terms of the solution to the magnetizing branches.The algorithm is useful in numerical simulation of the transients in the electric networks including ferromagnetic components.A better agreement is shown between calculated and experimental results, which verifies the validity of the algorithm.

Figure 3 :
Figure 3: An electric network including N magnetizing branches.

Figure 7 :
Figure 7: An electric network containing two groups of magnetizing branches.

Figure 8 :
Figure 8: Single phase network taking matrix form.

e t 1 .Figure 9 :
Figure 9: a Transient equivalent circuit of R-L coupled branches; b transient equivalent circuit of transmission lines.