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

Authors: Xiaoqing Zhang
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

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.


---

## Body

## 1. 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 linkage-magnetizing 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.


## 2. 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 resistanceRm [6–9], as shown in Figure 1. The constant resistance Rm 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(2.1)λ=Ai+Btan-1iW,
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)(2.2)Ld=dλdi=A+BWW2+i2.
In the unsaturation region, the small magnetizing current yields a large linear inductance(2.3)Ld0≈A+BW.
In the saturation region, the magnetizing current becomes large and the dynamic inductance is reduced to(2.4)Lds≈A.
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.Figure 1
Circuit representation of magnetizing branch.Figure 2
Saturation curve.


## 3.  Algorithm for Transient Calculation

An electric network includingN magnetizing branches is depicted in Figure 3, where the N core loss resistances (Rm) have been incorporated into the linear network. Let λ and i denote the flux linkage and magnetizing current vectors, respectively, that is, λ=[λ1,λ2,…,λN]T and i= [i1,i2,…,iN]T. In order to obtain the solution to the N nonlinear inductances, 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 u0(u0=[u10,u20,…,uN0]T) and input impedance matrix Z0(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 u0 and Z0, 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(3.1)dλdt=u0-Z0i.
The above derivative can be changed into(3.2)dλdt=dλdididt=Ld(i)didt.
The nonlinear dynamic inductance matrix Ld(i) is given as(3.3)Ld(i)=diag⁡[dλ1di1,dλ2di2,…,dλNdiN],
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 Ld(i) can be determined by (2.2). Substituting (3.2) into (3.1) leads to (3.4)didt=Ld-1(i)(u0-Z0i).
Letting F(i)=Ld-1(i)(u0-Z0i), this can be rewritten in a simplified form(3.5)didt=F(i).
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)in+1=in+k1H1+k2H2+⋯H1=Δt[F(in)+a1A(in)H1]H2=Δt[F(in+b1H1)+a2A(in+c1H1)H2]⋮⋮,
where A is the Jacobian matrix, A(i)=∂F(i)/∂i. Through a quite tedious calculation, it is possible to expand in+1-in 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)Δt:k1+k2-1=0(Δt)2:k1a1+k2(a2+b1)-12=0(Δt)3:{k1a12+k2[a22+(a1+a2)b1]-16=0k2(a2c1+12b12)-16=0.
These relationships involve six adjustable constants, k1, k2, a1, a2, b1,  and c1. Choosing the first two constants k1=0.75 and k2=0.25 [13], the remaining four constants can be determined as(3.8)a1=a2=12(1+13)=0.788675134595b1=-23=-1.15470053838c1=0.
Thus, (3.6) is simplified as a two-stage and third-order difference scheme(3.9)in+1=in+0.75H1+0.25H2,H1=Δt[1-0.788675134595ΔtA(in)]-1F(in),H2=Δt[1-0.788675134595ΔtA(in)]-1F(in-1.15470053838H1),
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 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].Figure 3
An electric network includingN magnetizing branches.Figure 4
Thevenin equivalent procedure.Figure 5
Thevenin equivalent network expressed in matrix form.Figure 6
Completely linear network.


## 4. Case Study

An electric network with two groups of magnetizing branches at both ends of transmission lines is shown in Figure7, 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: (4.1)e(t)=[1.0cos⁡(ωt+180∘),1.0cos⁡(ωt+60∘),1.0cos⁡(ωt-60∘)]Tp.uZc=diag⁡[606.87,279.88,279.88]Ω,τ=[1.43×10-3,1.01×10-3,1.01×10-3]Ts,l=305km,L=[0.592-0.149-0.149-0.1490.592-0.149-0.149-0.1490.592]H,R=[13.08-3.27-3.27-3.2713.08-3.27-3.27-3.2713.08]Ωλt=0.015it+0.691tan-1it0.142p.u.,Rmt=83.839kΩ,λr=0.31ir+0.83tan-1ir0.91p.u.
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, RRL and IRL(t) are the equivalent resistance matrix and current source vector of the R–L coupled branches, and Z, Ils(t) and Ilr(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.Figure 7
An electric network containing two groups of magnetizing branches.Figure 8
Single phase network taking matrix form.(a) Transient equivalent circuit ofR–L coupled branches; (b) transient equivalent circuit of transmission lines.
(a)(b)Figure 10
Equivalent calculation network.Calculated and experimental voltage waveforms.
(a)(b)(c)Contrast voltage waveforms.
(a)(b)(c)


## 5. 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.



---

*Source: 479626-2011-05-10.xml*