this page or responsible for its content.
Magnetic Cores Modeling for Ferroresonance Computations using the Harmonic Balance Method
1

Abstract-- The determination of the risk of ferroresonance in
actual HV or MV networks and the design of damping devices
need the use of accurate models. In this frame, the nonlinear
parts of the circuits, i.e. the magnetic cores, require a special
attention. This paper describes a modeling of the magnetization
curve and the core losses appropriate for ferroresonance
computations using the harmonic balance method. Tests
performed in order to get the parameters are discussed.

Index Terms Ferroresonance, Harmonic balance method,
High Voltage, Hysteresis, Magnetization, Modeling, Saturation,
Voltage Transformer

I. I
NTRODUCTION

ERRORESONANCE is a nonlinear phenomenon occur-
ring in electrical circuits involving at least one or several
saturable reactors, capacitances and a power supply.
Often, the saturable reactors are inductive potential
transformers. This phenomenon is characterized by the
possible existence of several stable regimes. These regimes
may be either periodic, with a base frequency equal to the
power supply frequency or to a sub multiple of it, pseudo-
periodic or chaotic.

The harmonic balance method, which is a particular case
of the Galerkin method, has proven powerful for the study of
periodic, but also pseudo-periodic regimes, in the above
mentioned circuits [e.g. 1-5]. For this method, a limited
Fourier series is used to represent the periodic (or pseudo-
periodic) behavior of the state variables. For instance, the
expression of the magnetic flux in the magnetizing branch of
the nonlinear inductance is :
( )
(
)
(
)
(
) + + =
K
k
s
k
c
k
t
k
t
k
t sin
cos
,
,
0
(1)
The corresponding Fourier coefficients for the current in this
branch are given by integration from the magnetic
characteristic
( ) i
of the inductance :

N. Janssens is with ELIA, Belgian Transmission System Operator, 125,
Rodestraat, 1630 Linkebeek, Belgium (e-mail:
noel.janssens@elia.be
) and
with the University of Louvain at Louvain La Neuve, Belgium.

( )
( ) (
)
dt
t
k
t
i
T
I
T
c
k
cos
2
0
, =
(2)
and with similar expressions for the sine terms
s
k
I
,
and for
the DC component
0
I
. The coefficients
c
k
I
,
,
s
k
I
,
,
0
I
are
nonlinear functions of all the coefficients
c
k
, ,
s
k
, ,
0 .
The harmonic balance method consists in introducing the
limited Fourier series in the differential equation of the circuit
and forcing to zero the contributions to each considered
harmonic component. So, an algebraic set of nonlinear
equations in the Fourier coefficients is obtained and may be
solved by using a general purpose routine.

The linear part of the circuit may be represented by its
Thevenin equivalents (voltage source
k
E
and complex
impedance
k
Z
) for the different frequencies k of the set
K
.
The use of Thevenin equivalents allows reducing the number
of equations to be solved : there is only one (complex)
equation for each harmonic component for each nonlinear
component, instead of one equation for each harmonic
component for each reactive element (linear or nonlinear). It
also facilitates the choice of an initial approximation of the
solution to be introduced in the computation program.

The use of the harmonic balance method may be extended
to find directly the domain limits in some parameter space of
the various kinds of ferroresonant regimes. This extension
consists in adding one equation to the above mentioned
algebraic set equating to zero the determinant of the Jacobian
of the harmonic balance equations relative to a set
K
(not
necessarily identical to the set
K
) [1-4].
.
II. B
ASIC MODELING PRINCIPLES

To get accurate results when studying a practical situation,
all the relevant components of the system must be modeled
adequately. In particular, the nonlinear reactors need a special
attention because the ferroresonance phenomenon pushes
these elements in high saturation beyond the validity domain
of the usually available data. The present paper is devoted to
the presentation of a saturable reactor modeling able to get
Magnetic Cores Modeling for Ferroresonance
Computations using the Harmonic Balance
Method

N. A. Janssens, Senior Member, IEEE
F
2
quantitative accurate results when used to study practical
situations.
The iron core magnetization characteristic has of course a
great influence, since the ferroresonance phenomena are
strongly related to this nonlinear function. Very often, the
question to be answered for estimating the risk of
ferroresonance or designing a damping circuit is : what is the
lower bound of the voltage source for which a given kind of
ferroresonant regime exists ? This question has an energetic
background : to what extent can the voltage source bring
enough energy to compensate the losses of a ferroresonant
regime ? From computation and tests, it may be concluded
that the system losses have a great importance on the lower
limit of the voltage source interval where a given
ferroresonant regime exists. A precise modeling of these
losses is of prime importance. It is to be noted that, when the
voltage is not close to these interval limits, the waveforms are
not very sensitive to these losses. Consequently, a validation
solely based on waveform comparisons is not very relevant.
For the time simulation of the system, a static magnetic
hysteresis model, like those described in [6,7] may be used,
associated with a conductance to represent the eddy current
losses and series resistances and leakage inductances. Such a
model could also be used for the harmonic balance method
and proved to provide good results for a low voltage
laboratory test [8]. However, the computation time may be
greatly reduced by using a univocal function
( ) i
to represent
the magnetic characteristic. In this case, the terms of the
Jacobian may be expressed by the harmonic components of
( )
( )
( )
( ) d
t
i
d
t
i
= [8]. For instance, let us consider the term
c
k
c
k
I
,
2
,
1 with
0
1 k
,
0
2 k
,
2
1
k
k . Using (1), we obtain
successively :
(
)
dt
t
k
d
i
d
T
I
c
k
T
c
k
c
k 1
cos
2
,
2
0
,
2
,
1 =
( ) (
) (
)
dt
t
k
t
k
i
T
T 1
cos
2
cos
2
0
=
c
k
k
c
k
k
I
I
,
2
1
,
2
1
2
1
2
1 + + =


Section IV will pay attention to the measurement of the
magnetic characteristic.
In order to have a clean convergence of the iteration
process, the function
( ) i
must be continuous. On the other
hand, the modeled magnetic characteristic must be very close
to the measured curve. For this purpose, an analytical
expression like, for instance, a polynomial will, in general,
show discrepancies with respect to the real curve or exhibit
oscillations in the (
i
) plane. Therefore, we found more
adequate to use a parabolic spline, consisting of successive
parabola segments with a continuous slope at the nodes.
Section IV will show an example of such a modeling.
Besides the magnetization curve, another important aspect
is the modeling of the core losses. During the computation
process to solve the algebraic system of equations, at the
beginning of each iteration, a set of Fourier coefficients is
given from the previous iterations (for the first iteration, the
user or an auxiliary routine provides an initial point). This
gives a flux (and voltage) behavior of the core magnetizing
branch. From there, it is possible to determine the value of
two linear conductances
F
G
and
H
G
such that the losses in
these conductances are equivalent respectively to the eddy
current losses and the hysteresis losses according to an
appropriate model. Doing so, the waveforms will differ
slightly from those obtained using a model to be used for a
time simulation. However, these differences will be very small
considering that the cores are made of soft magnetic materials
with limited losses and that the saturation has a much larger
effect on the waveform. The important point is that, globally
for the whole oscillation period, the core losses are accurately
modeled. The computation of the conductances
F
G
and
H
G

from basic data is developed in section III and an example of
core losses as a function of the saturation is shown in section
IV.
III. M
AGNETIC CORE LOSSES

Experimental studies have shown that the power
W

dissipated in silicon steel plates, for a distorted wave and in
the frequency range under interest here, can be written as the
sum of two terms relative to the hysteresis losses
H
W
and the
eddy current losses
F
W
:
F
H
W
W
W
+
=
(3)
of the form :
(
)
=
loops
H
H
w
f
W
max
min
,
(4)
(
)
rms
F
F
U
W
W
=
(5)
The symbols
W
(upper case letters) designate the powers
(energy per second) while the w (lower case) designate the
energy dissipated per cycle. In (4), the sum deals with all the
hysteresis loops (major and minor),
f
is the fundamental
frequency of the oscillation,
min and
max are the extreme
values of the flux. In (5),
rms
U
is the r.m.s. voltage across the
magnetizing branch.

A. Eddy current losses
We suppose that the function
0
F
W
giving the eddy
currents losses for a sine wave of the flux (directly rela