## Jelajahi eBook

## Kategori

## Jelajahi Buku audio

## Kategori

## Jelajahi Majalah

## Kategori

## Jelajahi Dokumen

## Kategori

0 penilaian0% menganggap dokumen ini bermanfaat (0 suara)

8 tayangan9 halamanIntegrasi dari histerisis untuk trafo asimetris dan dynamic modeling

Paper Referensi TA Histerisis pada Trafo Asimetris

© © All Rights Reserved

PDF, TXT atau baca online dari Scribd

Integrasi dari histerisis untuk trafo asimetris dan dynamic modeling

© All Rights Reserved

0 penilaian0% menganggap dokumen ini bermanfaat (0 suara)

8 tayangan9 halamanPaper Referensi TA Histerisis pada Trafo Asimetris

Integrasi dari histerisis untuk trafo asimetris dan dynamic modeling

© All Rights Reserved

Anda di halaman 1dari 9

Research Article

ISSN 1751-8660

Integration of the hysteresis in models of Received on 27th September 2015

Revised on 12th March 2016

asymmetric three-phase transformer: Accepted on 24th March 2016

doi: 10.1049/iet-epa.2015.0476

finite-element and dynamic electromagnetic www.ietdl.org

models

Faouzi Aboura 1,2 ✉, Omar Touhami 1

1

Laboratoire de Recherche en Électrotechnique, École Nationale Polytechnique, El Harrach 16200, Algeria

2

Energy Management Sector-Division, Siemens-Spa, 16035 Hydra, Algeria

✉ E-mail: faouzi.aboura@siemens.com

Abstract: This paper deals with the models elaboration for the asymmetric three-phase transformers. The first model uses

the principle of the hysteresis integration in two-dimensional finite-element method (FEM) with magnetic field

calculations. This method is based on inverse Jiles–Atherton hysteresis model. The second one is a dynamic

electromagnetic model (DEM) based on Tellinen hysteresis model. The DEM and FEM models used in ferromagnetic

cores have been modified by replacing the anhysteretic curve by the hysteresis loop. A simulation in time domain is

carried out and the simulated results are compared with those obtained experimentally to validate the proposed

approach. An error calculation is performed to show the accuracy of results obtained by DEM and FEM models.

ﬁeld parameter) [8].

The transformer is a static electrical machine that has been The hysteresis model used in the FEM keeps the same features that

represented by many models since more than a century. those of J–A hysteresis model original [9]. It has also the advantage

Nevertheless, the modelling accuracy of its operation has needed of being more suitable for the magnetic ﬁeld calculations with a

the integration of hysteresis phenomena in the models. For this magnetic vector potential formulation where magnetic induction is

reason, the modelling of hysteresis remains a subject of signiﬁcant beforehand known [10]. To solve the obtained non-linear

scientiﬁc interest [1, 2]. A great variety of hysteresis models are equations, we used the Newton–Raphson (NR) method leading to

available in the literature Jiles–Atherton (J–A), Preisach, Bouc– a differential reluctivity tensor. A DEM based on the duality

Wen model, Tellinen, Stoner–Wohlfarth, with their advantages and principle that considers the classical electrical equivalent circuit

their disadvantages. A comparison between some of these models coupled to a magnetic equivalent circuit has been implemented.

is given in [3]. The Tellinen hysteresis model that is a reversible scalar model was

In their representation, most models of transformers estimate the also used. The comparison between the simulated results with

saturation by an anhysteretic curve. This is mainly due to those obtained experimentally is presented in Section 5.

difﬁculties of correctly representing the hysteresis [2]. This led us The main contribution of this paper is the introduction of the

to introduce the hysteresis loops into the ferromagnetic cores of hysteresis in both models: FEM and DEM. The inverse J–A vector

the models of asymmetrical three-phase transformers. hysteresis model using an external user function was introduced in

This kind of transformer is called an asymmetric transformer due FEM software. The Tellinen model was also introduced in DEM.

to the difference in the length of the centre limb that is shorter than Unlike existent methods, this last model is developed in the

those of the other two limbs. Matlab/Simulink program and it includes all losses by hysteresis,

The aim of this paper is to provide an accurate dynamic and by eddy current in the lamination. It also includes a dynamic

electromagnetic model (DEM) for the representation of the hysteresis behaviour and not only quasi-static. The interest of this

asymmetric three-phase transformer. model lies in representation of all the parameters of the

For this purpose, we have organised our paper as follows. transformer, required to transient in low frequency. It can also be

Finite-element method (FEM) model and DEM based on duality coupled with FE model, which considers the current in the

principle and the dynamical hysteresis of Tellinen are presented in winding as an input quantity [11].

Sections 2 and 3. Then simulation in a time domain was

performed to compare the results of both models with

experimental ones. FEM is known as valuable tool for solving 2 Description of FE model of the asymmetric

various electromagnetic problems in electrical machines, three-phase transformer

nevertheless the integration of hysteresis in ﬁnite-element (FE)

computations remains challenging. One of major difﬁculties The formulation of classical magnetic vector potential of the FEM is

resides in the fact that the basic variable of the FE formulation adopted. The inverse J–A vector hysteresis model is chosen for the

does not coincide with the main variable of the hysteresis model. FEM. Particular attention is paid to the resolution of the non-linear

This problem may be solved by inverting the hysteresis model [4– equations by the NR method.

6]. The hysteresis in FEM model is based on the inverse J–A

hysteresis model. Its main problem is the perception that its ﬁve

parameters are all important. For example, if we consider α 2.1 Magnetic hysteresis in field equations

parameter when its value decreases, saturation, the remanent

induction and the coercive magnetic ﬁeld decrease [7]. It is well The hysteresis is a complex phenomenon related to physically

known that in many commercial software that use the J–A irreversible processes. Indeed, at a given moment a value of the

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

614 & The Institution of Engineering and Technology 2016

magnetic material properties depends not only on the intrinsic where [Vr] is the tensor of the relative reluctivity of the material; V0 is

properties of the material, but also on its history. There are an the reluctivity of the vacuum, and its expression is

inﬁnite number of possible curves B(H ), which makes it difﬁcult

the modelling of the hysteresis loop and consequently its 1 1

introduction into a FE program. V0 = = (4)

m0 4p.10−7

We have integrated the hysteresis in two-dimensional (2D)

time-domain solver FLUX [12], by using external user function A is the vector potential (Wb/m); Hc is the coercive magnetic ﬁeld;

named user subroutines. User subroutines are written in Groovy, [s] is the tensor of the conductivity (in S); V is the electric scalar

this is the name of an object-oriented programming language potential (in V).

intended to the Java platform that implements the inverse J–A The state variables are the magnetic vector potential A and the

vector hysteresis model. However, there exist several ways to electric scalar potential V [16].

implement the hysteresis model in ﬁeld calculation [13]. The differential reluctivity is deﬁned as

FE programs were traditionally developed in the Fortran and C

languages, which support procedural programming. During the last

dH DH H(t + Dt) − H(t)

15 years, the development of FE has gradually shifted towards an ∂v = = (5)

object-oriented approach. When the user subroutines are written in dB DB B(t + Dt) − B(t)

Fortran, we have to compile them using a Fortran compiler, and

when user subroutines are written in Groovy it uses Java syntax where dH and dB are the magnetic ﬁeld intensity and density,

and is directly compiled into bytecode by a Java compiler. The respectively.

Java is a simple language (simpler than C++) and efﬁciency of

Java code is comparable to that of C or C++ [14]. 2.2 Non-linear core model

Under rotational ﬂux observed in the corner joints (T-joints) of the

asymmetric three-phase transformer, the vector relationship between The J–A scalar hysteresis model is a model, in which the whole

magnetic ﬁeld and induction must be considered [15]. magnetisation is decomposed into reversible component Mrev

The differential reluctivity in the 2D modelling of an which corresponds to domain wall bending during the

asymmetrical three-phase transformer is introduced. The vector magnetisation process, and irreversible component Mir which

model is available by default in the FE software and we must take corresponds to domain displacement against the pinning effect [17]

that into account. This is why we used the inverse J–A model

instead of Tellinen model that is a scalar model. The main M = Mrev + Mir (6)

equations of FEM model are as follows.

Maxwell–Faraday equation

The model is written as a ﬁrst-order ordinary differential equation.

The main equations of J–A hysteresis model are presented in [18]

−∂B

rot(E) = (1)

∂t dM (1 − c) dMir /dBe + c/m0 dMan /dHe

=

dB 1 + m0 (1 − c)(1 − a) dMir /dBe + c(1 − a) dMan /dHe

implies the presence of the electric scalar potential V, such as

(7)

∂A dMrev = c(dMan − dMir ) (8)

E=− − gradV (2)

∂t

He = (H + aM) (9)

is

Be = m0 He (11)

∂A dMir Man− Mir

div(V0 [V r ]rot(A) − H c ) + [s] + gradV = 0 (3) = (12)

∂t dBe m0 k d

Fig. 1 DEM for wye/wye asymmetric three-phase transformer core with Tellinen hysteresis model

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

& The Institution of Engineering and Technology 2016 615

Fig. 2 Matlab/Simulink subsystem of asymmetric three-phase transformer

where μ0 is the vacuum magnetic permeability and Ms, a, c, α, k are M is the total magnetisation. k, c and a are second rank tensors which

the model parameters to be estimated, while δ is given by δ = sign terms must be obtained experimentally and j is the diagonal matrix

(dH/dt). He is the effective magnetic ﬁeld and Be is the effective of the anhysteretic functions derivatives against the effective ﬁeld

magnetic induction. Man is the anhysteretic function. This function component. Knowing dM we can write dH as follows:

can have different shapes depending on the type of each hysteretic dH = ∂v · dB

system and, until today, the Langevin function is the most often Where ∂v is the differential reluctivity tensor. For the 2D case,

used anhysteretic function in the J–A hysteresis model [19] the differential reluctivity tensor assumes the form of

H a ∂ Hx ∂Hx

Man = Ms coth e − (13)

a He ∂B ∂By ∂vxx ∂vxy

x

∂v = = (17)

∂ Hy ∂Hy ∂vyx ∂vyy

From a vector generalisation of the J–A scalar hysteresis model and

∂Bx ∂By

in its inverse version, the main equation of the model becomes [10,

20]

where x and y are the components of orthogonal axes.

1 All the tensor terms in (17) are given in detail in [10]. The

dM = [1 + fx (1 − a) + cj(1 − a) ]−1 ·[ fx + cj] dB (14) parameters of the J–A inverse method are the same as those of the

m0

−1

fx = xf xf xf (15)

Fig. 3 Matlab/Simulink subsystem for A, B, C in Fig. 2 Fig. 4 Family of the dynamic hysteresis loops

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

616 & The Institution of Engineering and Technology 2016

Table 1 Dynamical Tellinen hysteresis model parameters Table 2 Specification of asymmetric three-phase transformer 3 kVA, 50

Hz

Phase 1 Phase 2 Phase 3

Three-phase transformer parameters Value

α1 = 7.808 × 10−4 α2 = 7.808 × 10−4 α3 = 7.808 × 10−4

β1 = 0.3214 β2 = 0.3214 β3 = 0.3214 rated primary voltage, V 220

σ1 = 9.52 σ2 = 4.76 σ3 = 9.52 rated secondary voltage, V 127

σe(1) = 26.25 σe(2) = 13.125 σe(3) = 26.25 high-voltage windings turns 250

low-voltage windings turns 153

RP, Ω 0.6566

Rs, Ω 0.72466

a, A/m 22.05

α 9.22 × 10−6

k, A/m 10.62

Ms, A/m 1,810,000

c 0.15

the ‘branch-and-bound’ optimisation algorithm used in [17]. Fig. 6 Measurement circuit on non-linear (l–I) characteristics for limbs 1

and 2

3 Description of the DEM of the asymmetric losses in the transformer. The complete model of the

three-phase transformer electromagnetic circuit is shown in Fig. 1.

The detailed description of the subsystem ‘asymmetric

The model uses electric and magnetic coupling, the non-linear core three-phase transformer’ is shown in Fig. 2. Flux linkages and

representation, together with the dynamical behaviour due to ﬂuxes are not mixed and the leakage air ﬂux path is modelled by a

integration of eddy current factor in the ﬁeld equations. The single linear reluctance term obtained from the zero-sequence tests

hysteresis model is described in detail in [21]. It is ﬂexible [22]. The connection with the dynamical formulation is given by

because we can use the variation of the magnetic induction to the subsystems A, B, C as shown in Fig. 3, through a Matlab

compute the magnetic ﬁeld and vice versa. This model can function.

simulate the quasi-static hysteresis loops, as well as the dynamic The Matlab function describes (18)–(23) and the global hysteresis

hysteresis loops, that consider the eddy currents and the excess behaviour. The controlled MMF source (1, 2 and 3) across the

Fig. 7 Experimental setup for three-phase transformer tests to measure the magnetic characteristic (l–I)

a Side view of the experimental setup

b Front view of the experimental setup

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

& The Institution of Engineering and Technology 2016 617

subsystem (1, 2 and 3) is used to compute (21) and (22) based on

(18)–(20). At every time step (Δt = 1 × 10−6), an ‘IF statement’

selects one of (21) and (22) based on the sign of the derivative of

ﬂux. The script of the Matlab function of the subsystems 1 and 3

is given in Appendix.

To take into account of core asymmetry of the three-phase

transformer type, each core limb is modelled with its own

customisable non-linear hysteresis. The main equations of the

Tellinen model use the (l–I ) characteristics of all three limbs that

are different. They are described as

+

dwx (fx ) a ·b

= x x (18)

dfx bx fx + sx + 1

w+

x (fx ) = sgn(fx − sx ) · sx loge (bx fx − sx + 1) (19)

w−

x (fx ) = sgn(fx + sx ) · sx loge (bx fx + sx + 1) (20)

w+ −

x , wx are, respectively, the ascending and descending parts of the

limiting hysteresis curve. jx is the instantaneous limb ﬂuxes, for

each phase x, dependent on the limb MMF potential fx.

ρx is the slope of the fully saturated region along the limiting

hysteresis curve. The hysteresis equations (21) and (22) are given by

dwx dfx w− (f ) − wx (fx ) dw+ x (fx )

= rx + −x x · − r (21)

dt dt wx (fx ) − w+

x (fx ) dfx x

If (djx/dt) ≥ 0

−

dwx dfx wx (fx ) − w+

x (fx ) dwx (fx )

= r + · − rx (22)

dt dt x w− +

x (fx ) − wx (fx ) dfx

This function has the advantage of requiring only two ﬁtting

parameters (αx, βx) to shape the curve to measured experimental

loops. αx and βx inﬂuences the vertical and horizontal scaling,

respectively, of jx( fx). σx controls the width of the major

hysteresis loop. The measured hysteresis loop of each limb (Fig. 8)

is used to measure the two ﬁtting parameters (αx, βx), and σx

controls the width of the major hysteresis loop for each limb (x =

1, 2, 3). In addition to (18)–(22), we used the eddy current factor

σe(x) for the lamination in (23). This factor reﬂects the classical

eddy current and the excess losses and fxd is the applied dynamic

current [21]

dwx

fxd = fx + se(x) (23)

dt

shown in Fig. 4. The Tellinen dynamic hysteresis model

parameters are presented in Table 1.

When σe(x) is equal to zero, then hysteresis is quasi-static as shown

in Fig. 5. This family of quasi-static hysteresis loops has the same

Fig. 8 Hysteresis loop of the limbs 1, 2 and 3 obtained from

magnetic coercive ﬁeld.

experimentations

a Hysteresis loop of the limb 1

b Hysteresis loop of the limb 2

4 Experimental measurement of hysteresis loops c Hysteresis loop of the limb 3

with rated power at 3 kVA, and a voltage primary and a secondary Nevertheless, we have used the method of Fuchs [22] for

are, respectively, rated at 220 and 127 V, operated at 50 Hz. The measuring the hysteresis curve in each limb of the three-phase

transformer is supplied from a variable three-phase power source transformer, and the cross-coupling between different limbs due to

using an autotransformer. The data of the asymmetric three-phase saturation is taking place in the magnetic circuit of the DEM.

transformer are presented in Table 2. To measure the magnetic characteristic (l–I) for each limb of the

The currents and voltages are recorded through data acquisition transformer it is important to excite two phases of the three-phase

devices via current and voltage sensors. winding and to have two limbs to have the same ﬂux magnitudes

A three-phase core-type transformer iron core model by taking in opposite directions (as shown in Fig. 6), i.e. two phases of the

into account the effects of cross-couplings between different limbs low-voltage windings are connected in parallel with reverse

due to saturation has been published by Dolinar et al. [23]. polarity. This method was deﬁnitively conﬁrmed in another paper

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

618 & The Institution of Engineering and Technology 2016

Fig. 9 M (H) characteristic of the core material ET130-35

‘measurement of (l–I ) characteristics of asymmetric three-phase

transformers and their applications’ [24]. Fluxes (l1) and (l2) can

be calculated by integrating the induced voltages (V1) and (V2)of

the secondary windings, respectively, of the two (1, 2) excited

phases by an analogical RC circuit. The laboratory test transformer

to measure the magnetic characteristic (l–I ) for each limb of the

transformer is shown in Fig. 7.

The magnetic characteristics (l–I ) of all limbs are different

(Fig. 8) because they depend not only on magnetic properties of

core material but also from geometric dimensions of each limb,

consisting of vertical part and horizontal portion.

those obtained experimentally

and FEM models, with those obtained from experimental tests.

FEM model. Some transformer characteristics are listed as follows

(primary/secondary): 220/127 V, 250/153 turns. The hysteresis

characteristics of the material are shown in Fig. 9.

The transformer core is made of 0.35 mm Fe-3%Si sheets. The

non-linear FEM equations are solved by means of the NR method.

A relative tolerance of 10−5 is adopted.

The 2D automatic mesh with triangular elements on the

three-phase transformer and the study domain are shown in Fig. 10.

Fig. 11 Comparison between experimental no-load currents in steady state

and the FEM ones

a Steady-state current in phase 1

b Steady-state current in phase 2

c Steady-state current in phase 3

steady state with those obtained by the simulated FEM model

with the inverse J–A vector hysteresis model is shown in

Fig. 11. The FE model due to high number of parameters and

its accuracy in calculations requires a very important time of

calculation. In addition, the determination of the parameters of

the J–A model needs precise knowledge of the magnetic

material of the transformer. All the results are as observed on

our notebook with the following conﬁguration: Windows 7,

CPU: Intel(R) Core(TM) i5-2410M @ 2.3 GHz and RAM:

6 GB and the time duration to obtain 1 period of current

Fig. 10 Two-dimensional FE meshes with triangular elements signal (20 ms) is 36 min.

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

& The Institution of Engineering and Technology 2016 619

Fig. 12 Comparison between experimental no-load currents in steady state

and the DEM ones Fig. 13 Comparison between the inrush currents in DEM and the

a Steady-state current in phase 1 experimental ones

b Steady-state current in phase 2 a Inrush current in phase 1

c Steady-state current in phase 3 b Inrush current in phase 2

c Inrush current in phase 3

determined precisely. The comparison between the inrush currents of

The analysis with DEM does not require an important calculation DEM and the experimental ones is presented in Fig. 13. The inrush

time. It features a large convergence and stability. It is also being current of a transformer depends also on the residual ﬂux of the

adapted to calculations of the electromagnetic transients, such as: magnetic core [25]. The slightest phase mismatch can lead to large

the inrush current and the ferroresonance. discrepancies.

The comparison of the experimental no-load currents in steady

state with those simulated by DEM is shown in Fig. 12. 5.4 Ferroresonance study

5.3 Inrush currents We want clearly show the validity of the developed model through a

ferroresonance test. The realised circuit is described in [26]. We

To obtain good comparison of measured and simulated waveforms, chose the fundamental mode because of the ease of obtaining it

the initial conditions (i.e. the instant of voltage phase angle) must be compared with other modes. The obtained currents can be used as

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

620 & The Institution of Engineering and Technology 2016

† The developed DEM is well suitable for the study of the

low-frequency transients. In addition, the computation time in the

DEM is very short.

7 References

1 de Léon, F., Farazmand, A., Joseph, P.: ‘Comparing the T and II equivalent circuits

for the calculation of transformer inrush currents’, IEEE Trans. Power Deliv., 2012,

27, (4), pp. 2390–2398

2 Moses, P.S., Massoum, M.A.S., Toliyat, H.A.: ‘Dynamic modelling of three-phase

asymmetric power transformers with magnetic hysteresis: no-load and inrush

currents’, IEEE Trans. Energy Convers., 2010, 25, (4), pp. 1040–1047

3 De Leon, F., Semlyen, A.: ‘A simple representation of dynamic hysteresis in power

transformers’, IEEE Trans. Power Deliv., 1998, 10, (1), pp. 315–321

4 Gyselinck, J., Vandevelde, L., Melkebeek, J.: ‘Calculation of no load losses in an

induction motor using an inverse vector Preisach model and an eddy current loss

model’, IEEE Trans. Magn., 2000, 36, (4), pp. 856–860

5 Sadowski, N., Batistela, N.J., Bastos, J.P.A., et al.: ‘An inverse Jiles–Atherton

model to take into account hysteresis in time-stepping ﬁnite-element

calculations’, IEEE Trans. Magn., 2002, 32, (2), pp. 797–800

Fig. 14 Comparison of the current I1 in fundamental mode between DEM 6 Kefalas, T., Kladas, A.: ‘Harmonic impact on distribution transformer no-load

and the experimental results loss’, IEEE Trans. Ind. Electron., 2010, 57, (1), pp. 193–200

7 Bui, A.T.: ‘Caractérisation et modélisation du comportement des matériaux

magnétiques doux sous contrainte thermique’. PhD thesis, Claude Bernard

University, 2011

Table 3 Calculation error for currents obtained by DEM and FEM 8 Marinucci, T., Maglietta, V.: ‘Identiﬁcation of the Jiles–Atherton parameters using

models with those obtained experimentally commercial software’, IEEE Trans. Power Deliv., 2004, 19, (1), pp. 200–207

9 Jiles, D., Atherton, D.: ‘Ferromagnetic hysteresis’, IEEE Trans. Magn., 1983, 19,

RMS error DEM FEM (5), pp. 2183–2185

10 Leite, J.V., Sadowski, N., Kuo-Peng, P., et al.: ‘Inverse Jiles–Atherton vector

I1, % 2.3096 3.3174 hysteresis model’, IEEE Trans. Magn., 2004, 40, (4), pp. 1769–1775

I2, % 2.5842 5.4670 11 Charalambous, C.A., Wang, Z.D., Jarman, et al.: ‘Two dimensional ﬁnite element

I3, % 3.2298 4.9083 electromagnetic analysis of an auto-transformer experiencing ferroresonance’,

IEEE Trans. Power Deliv., 2009, 24, (3), pp. 1275–1283

12 ‘Flux software user guide’, http://www.cedrat.com, accessed 27 January 2014

13 Mathekg, M.E., McMahon, R.A., Knight, A.M.: ‘Application of the ﬁxed point

inputs to the FEM and so show the inﬂuence of other ferroresonant method for solution in time stepping ﬁnite element analysis using the inverse

vector Jiles–Atherton model’, IEEE Trans. Magn., 2011, 47, (10), pp. 3048–3051

modes on the transformer [11, 27, 28]. The comparison of currents in 14 Nikishkov, G.: ‘Programming ﬁnite element in java’ (Springer-Verlag, London,

the fundamental ferroresonant mode of DEM with those 2010)

experimental is shown in Fig. 14. 15 Ferreira da Luz, M.V., Leite, J.V., Benabou, A., et al.: ‘Three-phase transformer

modeling using a vector hysteresis model and including the eddy current and the

anomalous losses’, IEEE Trans Magn., 2010, 46, (8), pp. 3201–3204

5.5 Comparison of no-load currents in steady-state 16 Bastos, J.P.A., Sadowski, N.: ‘Electromagnetic modelling by ﬁnite element

operation between those of the models and those of methods’ (Marcel Dekker, New York, 2003)

17 Chwastek, K., Szczygłowski, J.: ‘Modelling dynamic hysteresis loops in steel

experimental tests sheets’, COMPEL, Int. J. Comput. Math. Electr. Electron. Eng., 2009, 28, (3),

pp. 603–612

For the comparison of the obtained results from DEM and FEM 18 Lacerda Ribas, J.C., Lourenco, E.M., Leite, J.V., et al.: ‘Modelling ferroresonance

models with those experimental, we have used the criteria of the phenomena with a ﬂux-current Jiles–Atherton hysteresis approach’, IEEE Trans.

Magn., 2013, 49, pp. 1797–1800

root mean square (RMS) error, deﬁned by 19 Kokornaczyk, E., Gutowski, M.W.: ‘Anhysteretic functions for the Jiles–Atherton

model’, IEEE Trans. Magn., 2015, 51, (2), pp. 1–5

n 2 n 2 20 Leite, J.V., Benabou, A., Sadowski, N.: ‘Transformer inrush currents taking into

1 xref − 1 xsim

RMS =

n 2

(24) account vector hysteresis’, IEEE Trans. Magn., 2010, 46, (8), pp. 3237–3240

21 Tellinen, J.: ‘A simple scalar model for magnetic hysteresis’, IEEE Trans. Magn.,

1 xref

1998, 34, (4), pp. 2200–2206

22 Fuchs, E.F., You, Y.: ‘Measurement of l–i characteristics of asymmetric

where xref is experimental result and xsim is a simulated result. three-phase transformers and their applications’, IEEE Trans. Power Deliv.,

The comparison of no-load currents in steady state, simulated by 2002, 17, (4), pp. 983–990

FEM and DEM with those measured, shows that DEM gives better 23 Dolinar, M., Dolinar, D., Stumberger, G., et al.: ‘A three-phase core-type

transformer iron core model with included magnetic cross saturation’, IEEE

results (see Table 3). Trans. Magn., 2006, 42, (10), pp. 2849–2851

24 Fuchs, E.F., You, Y.: ‘Closure on measurement of – characteristics of asymmetric

three-phase transformers and their applications’, IEEE Trans. Power Deliv., 2003,

6 Conclusions 18, (2), pp. 642–643

25 Naghizadeh, R.A., Vahidi, B., Hosseinian, S.H.: ‘Modelling of inrush current in

transformers using inverse Jiles–Atherton hysteresis model with a neuro-shufﬂed

A methodology to include the magnetic hysteresis models in the FE frog-leaping algorithm approach’, IET Electr. Power Appl., 2012, 6, (9), pp. 727–734

calculations and DEM of an asymmetric three-phase transformer is 26 Moses, P.S., Masoum, M.A.S., Toliyat, H.A.: ‘Impacts of hysteresis and magnetic

presented in this paper. couplings on the stability domain of ferroresonance in asymmetric three-phase

three-leg transformers’, IEEE Trans. Energy Convers., 2011, 26, (2), pp. 581–592

The main conclusions we have reached are: 27 Charalambous, C.A., Wang, Z.D., Jarman, P., et al.: ‘Time-domain ﬁnite-element

technique for quantifying the effect of sustained ferroresonance on power

† A newly DEM based on Tellinen hysteresis model is developed. transformer core bolts’, IET Electr. Power Appl., 2014, 8, (6), pp. 221–231

28 Charalambous, C.A., Wang, Z.D., Jarman, P., et al.: ‘Frequency domain analysis of

The obtained dynamic behaviour from an eddy current factor is a power transformer experiencing sustained ferroresonance’, IET Gener. Transm.

included in a Matlab/Simulink program. Distrib., 2011, 5, (6), pp. 640–649

† The DEM and FEM models are numerically stable. They give

good performance in terms of numerical convergence.

† The results obtained on an asymmetric three-phase transformer 8 Appendix

are very satisfactory and conﬁrm the validity of the developed

models. 8.1 Implementation of subsystems 1, 3 in Matlab script

† The computation time in the FEM model with the external language

hysteresis function is very high by comparison to the case when

an anhysteretic curve is used. function [Imd,Im] = fcn(im,qm,dqm)

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

& The Institution of Engineering and Technology 2016 621

%#codegen dq_m = a*b./(b*abs(im + seg)+1);

mu = 4*pi*10^-7; if dqm > 0

s = 32e-4;l = 0.56; dim = dqm*1/abs((ρ + (q_m-qm)/(q_m-q_p)*(dq_p-ρ)));

ρ = mu × s/l; else

a = 0.0007808;b = 0.3214;seg = 9.52;segma = 26.25; dim = dqm*1/abs((ρ + (qm-q_p)/(q_m-q_p)*(dq_m-ρ)));

q_p = (im-seg)./abs(im-seg)*a.*log(b*abs(im-seg) + 1); end

q_m = (im + seg)./abs(im + seg)*a.*log(b*abs(im + seg)+1); Im = im + dim*1e-6;

dq_p = a*b./(b*abs(im-seg)+1); Imd = Im + segma*dqm;

IET Electr. Power Appl., 2016, Vol. 10, Iss. 7, pp. 614–622

622 & The Institution of Engineering and Technology 2016