A comprehensive finite-element model of a turbine-generator infinite-busbar system

Available online at www.sciencedirect.com

Finite Elements in Analysis and Design 40 (2004) 485 – 509 www.elsevier.com/locate/ nel

A comprehensive nite-element model of a turbine-generator in nite-busbar system Rafael Escarela-Pereza;∗ , Marco A. Arjona-Lopezb , Enrique Melgoza-Vazquezc , Eduardo Campero-Littlewooda , Carlos Aviles-Cruzd a

Departamento de Energia, Universidad Autonoma Metropolitana-Azcapotzalco, Av. San Pablo 180, Col. Reynosa, Azcapotzalco, Mexico, D.F. C.P. 02200, Mexico b Departamento de Ingenieria Electrica, Instituto Tecnologico de la Laguna, 27000 Torreon, Coah. Mexico c Departamento de Ingenieria Electrica, Instituto Tecnologico de Morelia, Mich. Mexico d Departamento de Electronica, Universidad Autonoma Metropolitana, Azcapotzalco, Mexico Received 1 September 2002; received in revised form 9 December 2002; accepted 23 December 2002

Abstract This paper shows the development of a two-dimensional nite-element magnetic model of a turbine generator coupled to an in nite busbar through a transmission line and a delta-star connected transformer. The nite-element equations of the machine are strongly coupled to the circuit equations of the transmission line and transformer. The circuit equations developed for the delta-star connection enable the simulation of several short circuit conditions at the transformer terminals with the aid of fault matrices. The calculation of electromagnetic torque is performed by surface integration of the Lorentz force, which leads to a unique result for the nite-element model used in this work. The coupling with the external controls of the machine (automatic voltage regulator and governor) is performed in a weak form. Nevertheless, any control model can be easily connected to the nite-element model. Since the use of this nite element model is meant for transient simulations, it is desirable to obtain some simpli cations to speed up calculations. Hence, a three-phase current sheet is proposed in this work to avoid the explicit modelling of the stator coreback. This approach also leads to a simple approach to model motion, where remeshing is not necessary. ? 2003 Elsevier B.V. All rights reserved. Keywords: Synchronous generators; Finite elements; Power system elements

∗

Corresponding author. Tel.: +52-55-5318-9584; fax: +52-55-5394-7378. E-mail address: [email protected] (R. Escarela-Perez).

0168-874X/$ - see front matter ? 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0168-874X(03)00074-X

486

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

1. Notation A symbol between square brackets denotes a matrix (e.g., [S]) while a symbol between braces denotes a vector (e.g., {A}). Subscript e refers to a single nite element. A ‘bar’ on top of symbols denotes a space vector (e.g., VF ). The following set of symbols will be used for the electromagnetic nite-element model: E H D B J A    V Az Jz I If J LeH P Jt Vf Pe Pm Vs !; !s 

Electric eld intensity Magnetic eld intensity Electric Gux density Magnetic Gux density Electric current density Magnetic vector potential Electric charge density Magnetic reluctivity Electric conductivity Electric potential diHerence z component of the magnetic vector potential z component of the electric current density Winding current Equivalent eld winding current at time t Rotor inertia Axial eHective length of the two-dimensional nite-element model Number of machine poles Time step Actual voltage supplied to the eld winding Electromagnetic air-gap power Mechanical power supplied by the prime mover In nite busbar voltage Rotor speed and synchronous speed Rotor angle with respect to the synchronous reference axis de ned by the phasor representation of Vs

2. Introduction The simulation of synchronous machines is an important subject of study in power system analyses and in the development of new techniques for their design. The rst attempts to describe the electromechanical behavior of these kind of machines can be traced back to the ending of the nineteen century. The representation of synchronous machines was greatly simpli ed with the development of the two-axis theory [1]. This allowed the practical calculation of synchronous machine performance without using digital computers (unavailable at that time). Power system prediction was greatly enhanced with the availability of the rst computers and several tests were devised to obtain the parameters of very simple two-axis equivalent circuits. However, it took several decades to note that

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

487

these circuits needed to be improved, since the eld winding circuit was poorly represented. The pioneering work of Canay [2] helped to obtain a much better approximation of the eld winding with the introduction of the leakage diHerential inductance concept. The two-axis equivalent circuits, used in modern power system calculations, have basically the same structure of Canay’s equivalent circuit, but more branches are employed to represent the rotor circuits [3]. The main problem associated with the two-axis approach is that the determination of parameters of high-order equivalent circuits is very diLcult since the tests that are designed for this purpose do not usually represent the actual operating conditions of the machine. The standstill frequency response test [4,5] is one of the most accepted tests for the determination of high-order circuit parameters due to its low risk of damage and reduced costs. Although it is now possible to obtain parameters of circuits of any order (for both the d and q axes) [6,7], there are several problems associated with local minima [8,9] and unrepresentative operating conditions [5]. Nevertheless, high-order equivalents circuits have opened up a much better representation of the synchronous machine. Another approach to model the synchronous machine is the Finite Element Method [10], which was rst applied by Chari [11] for the simulation of d-axis steady-state conditions. The method was then extended to incorporate load operating conditions, using additional constrains or incorporating phasor diagrams [12–15]. One of the rst solutions of synchronous machine time-dependent magnetic problems was given by Hannalla and Macdonald using the nodal method for the establishment of the nite-element equations [16]. These authors were also able to solve problems with unknown winding currents using conductivity matrices [17] for the rst time. The nite-element approach has since been complemented with diHerent techniques to allow for the connection of external circuits and mechanical systems [18–20]. The main advantage of nite-element techniques over traditional equivalent circuits is the rigorous incorporation of the machine geometry and the better modelling of nonlinear materials, leading to a very good match between test and simulation results. However, large computer resources may be necessary in transient simulations. This paper shows the full development of a nite-element model of a turbine-generator connected to an in nite busbar through a transmission line and delta-star connected transformer. Three current sheets are proposed to model the stator phases, in order to avoid direct modelling of the stator coreback, which is assumed to have an in nite permeability. This leads to a reduced nite-element mesh where more elements can be used to model the complex electromagnetic behavior of the solid rotors of turbine generators. The magnetomotive force drop of the stator coreback is taken into consideration by increasing the airgap length, using a simple but accurate approach based on open-circuit calculations. In addition, rotor motion is easily implemented with the stator current sheets, eliminating the need of remeshing and all of its associated numerical problems [21]. The rotor motion simulation only requires redistribution of stator currents, leading to a nite-element model which demands low computational resources and it is also very accurate. Incorporation of external electrical devices (connected to the machine windings) is performed using the conductivity [16,22] and winding vector [23] concepts, giving a strong coupling between the nite-element and circuit equations. Torque calculation is eLciently computed by surface integration of the Lorentz force within the stator conductors, as currents sheets lend themselves to obtain the torque value naturally with this technique. The connection of machine controls is performed in a weak form, but they can be easily included within the nite-element model. Hence, the nite-element model developed in this work is eLcient and easy to implement.

488

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

3. Finite-element model of synchronous-generator in!nite-busbar system 3.1. The 2D electromagnetic equations of the synchronous machine The solution of Maxwell’s equations [24] fully determines the electromagnetic behavior of synchronous generators and they are given in diHerential form by: ∇×E=−

9B ; 9t

∇×H=J+

9D ; 9t

∇ · D = ; ∇ · B = 0:

(1)

These equations are complemented with constitutive relationships [24] to take into consideration the material properties of the electromagnetic media involved in the machine. Although the numerical solution of these equations is possible, it is convenient to neglect some electromagnetic eHects in low-frequency power devices. One important simpli cation is obtained when the displacement current (9D=9t) is disregarded. As a result, it is assumed that power transfer from stator to rotor takes place instantaneously and that J represents the current densities impressed from external sources and the generated eddy currents on massive conductors. A further simpli cation is attained if a two-dimensional representation of the machine is accurate enough. This is particularly possible in large turbine generators since their axial lengths are much bigger than their diameters. Consequently, the machine can be reasonably considered as if it were of in nite axial length, as the magnetic eld distribution is nearly planar. Some three-dimensional eHects can be taken into consideration outside the nite-element model with external circuit/lumped parameters and others by proper adjustment of material properties (see below). The solution of the two-dimensional version of Maxwell’s equations is better obtained via the magnetic vector potential A [24], which has one component in the axial direction in two dimensions and is de ned by B = ∇ × A:

(2)

In order to completely de ne the magnetic vector potential, its divergence must be speci ed, which in two dimensions is automatically enforced to zero. By proper vector manipulation of Eqs. (1) and (2), it can be shown that the transient electromagnetic problem of the synchronous machine can be formulated with the following nonlinear diHusion equation: 
9Az 9V (∇ · ∇)Az =  + = −Jz ; (3) 9z 9t where the term (9V=9z) is the potential diHerence (V1 − V2 ) [20] measured at the far and near ends of conducting regions divided by the eHective length (LeH ) of the conductor. It is this potential diHerence that allows the coupling of the eld equations with the external devices connected

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

489

to the machine terminals, as well as the interconnection between diHerent internal regions of the electromagnetic model. Eq. (3), subject to proper boundary conditions, represents the solution of the transient two-dimensional electromagnetic problem in synchronous generators. Eq. (3) is more easily tackled by separating the eld model into three regions of diHerent conducting characteristics: (i) Current Free Regions. Here, Jz = 0, so Eq. (3) is simpli ed to Laplace’s equation: (∇ · ∇)Az = 0:

(4)

The airgap and the laminated ferromagnetic materials belong to this category. (ii) Massive Conductors, such as the solid steel forging, wedges and damper bars of the machine rotor. Eddy currents can exist in these regions. Eq. (3) is now given by: 9Az (5) 9t if it is assumed that these conducting regions are short-circuited at their far ends. The assumption of short-circuited regions at in nity is only an approximation, since rotor eddy currents have to pass through the rotor end bell (retaining ring) which has a nite impedance. Two ways are available to take into consideration this eHect. The rst method consists of coupling the rotor conductors with the end bell impedance by considering the (9V=9z) term in (5). In the second method, the rotor conductivities are properly increased to include the end bell eHect [18]. This latter method is chosen here since it eases the numerical implementation of the nite-element model. Nevertheless, it is important to emphasize that any external circuit can be coupled to the massive conductors of the machine model. (iii) Stranded conductors (or machine windings) are regions where the currents densities are assumed to be uniformly distributed over the area S occupied by a winding of N turns, such that the current density can be simply calculated as Jz = NI=S. This means that parasitic currents are negligible in these regions and that Eq. (3) can be rewritten as (∇ · ∇)Az = 

NI S for each machine winding. (∇ · ∇)Az = −

(6)

3.2. Finite-element versions of the nonlinear di8usion equation Although it would be preferable to obtain an analytical solution for the diHusion equation, this is not possible in most cases due to the complex geometries and nonlinear magnetic materials involved in synchronous generators. Nevertheless, very accurate solutions can be obtained with numerical techniques, such as the Finite Element Method [10]. Finite-element techniques consist of the discretization of the problem region into small subregions ( nite elements) where the unknown variable is approximated by a preselected function (polynomials are frequently chosen). The determination of values of the unknown function at discrete points (subject to known constraints) represents the solution of the boundary value problem. This discretization can be carried out using elements of any shape and order. Although higher order elements [10] may give more accurate results, rst-order triangular elements have been used in the analysis of low frequency electromagnetic devices for many years because they are easy to handle mathematically.

490

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

Furthermore, the intricate geometry of turbine-generators requires a great number of elements for the solution to be accurate. Hence, rst-order nite elements are chosen to construct the numerical model of the synchronous generator considered in this work. The numerical version of Eq. (4) in a single nite element is given by e [] {A}e = {0}e ; (7) 4Se e where {A}e is a 3 × 1 vector containing the potential values at element vertices and {0}e is a null vector. e and Se are the reluctivity and area of the element respectively. []e is a 3 × 3 matrix given by   b21 + c12 b 1 b2 + c 1 c2 b1 b 3 + c 1 c3   []e =  (8) b22 + c22 b2 b 3 + c 2 c3    b1 b2 + c 1 c2 b1 b3 + c 1 c3 b 2 b3 + c 2 c3 b23 + c32 whose coeLcients depend on the element node coordinates (xi ; yi ): a 1 = x2 y3 − x 3 y2 ; b1 = y2 − y3 ; c1 = x 3 − x 2 ;

(9)

a2 ; b2 ; : : : ; etc. are simply de ned by permutation of the subscripts. The nite-element version of Eq. (5) for massive conductors is found to be:

e 9A []e {A}e = −[J]e ; (10) 4Se 9t e where {9Az =9t}e gives the time derivative of the element nodal potentials and [J]e is a diagonal matrix with its entries calculated as ii = 13 Se e

i = 1; 2; 3;

where e is the element conductivity (constant for for stranded conductor regions yields: e [] {A}e = {I }e : 4Se e

(11) rst order elements). Discretization of Eq. (6) (12)

Here, {I }e is a vector of impressed nodal currents and its entries are calculated as a fraction of the total magnetomotive force NI existing in the stranded region of area S: Si N I i = 1; 2; 3; (13) S where Si is an area portion of the element that contains node i (usually one third of the element area). Once the elemental matrices have been established it is possible to merge them into a global matrix, modelling the whole electromagnetic problem region of the turbine-generator. Thus, Ii =

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

Eqs. (7), (10) and (12) are combined to give:

9A ; [S]{A} = [N ]{Iw } − [J] 9t

491

(14)

where {A} is a global vector of magnetic potentials at nodes, [S] embodies the element geometries and magnetic properties of the region being studied. This matrix is invariably symmetric and sparse. The matrix [J] (better known as conductivity matrix [16]) is formed by the material conductivities of the massive conductors. [N ] is a winding matrix that contains individual vectors known as winding vectors [23,25]. They allocate fractional numbers of winding turns to nite-element nodes that represent winding conductors (see Eq. (13)), leading to proper distribution of winding currents through the simple product of the winding vector with the winding current. Winding vectors have zero entries everywhere but in locations of element nodes representing the winding region. The sum of the entries of a winding vector always gives the total number of turns of the machine winding. Four winding vectors are de ned for a three-phase synchronous generator: [N ] = [{Na }{Nb }{Nc }{Nf }]

(15)

three for the stator phases a; b and c and one for the eld winding f. Winding vectors not only give the distribution of stranded conductor currents in the nite-element model but they can be also used to obtain the magnetic Gux linkages with machine windings. It is important to notice that the fraction of winding turns attributed to a particular node times the magnetic vector potential value at the node location gives the Gux magnetic linkage per unit length of the winding turns associated with that node. As a result, the product {N }T {A} gives the total Gux linkage, per unit length, with the winding. It is assumed here that only one pole pitch of the P pole machine is represented by the nite-element model. {Iw } is a vector that contains the four machine winding currents: {Iw } = {ia ib ic if }T :

(16)

Eq. (14) cannot be solved independently, if the problem is not current driven (a situation found in most practical cases). As a result, it is necessary to include the external device equations within the nite element formulation to have a consistent system of equations where the winding currents are unknown. 3.3. Coupling with external devices The synchronous machine model must be coupled to the external electrical devices connected to the machine: the excitation system and the electric power network. The circuit- eld coupled problem is solved here by using the idea that winding terminals within the nite-element model can be taken out and connected to the terminals of the external devices. This approach is numerically implemented through the nite-element versions of winding potential diHerences, leading to conductivity-like matrices. Alternative coupling methods based on similar principles are possible [19,20,26,27]. Nevertheless, they lead to diHerent structures of systems of equations and, therefore, diHerent methods of solutions may be necessary.

492

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

Fig. 1. Coupling of the nite-element winding terminals with the external excitation system.

3.3.1. Coupling with the excitation system Fig. 1 shows schematically the electrical connection between the nite-element model and the eld exciter circuit. Here, Lf and Rohf are the lumped parameters that approximate the three-dimensional/ overhang eHects of the eld winding. The potential diHerence given by the nite element model at terminals of machine windings can be calculated by adding the individual potential diHerences of each series connected conductor belonging to a particular winding. Then, if it is assumed that the conductivities of winding conductors are uniform over their cross sections, the entire winding potential diHerence for a two-pole machine can be evaluated as

9A T ek = −Ik Rk − 2LeH {Nk } ; (17) 9t where {Nk } is the winding vector of the kth winding. Machine windings constructed with parallel circuits can also be included with this formulation by replacing the original machine winding by an equivalent one that contains the same number of circuits but connected in series. Thus, each new series circuit has a number of turns which is determined by dividing the actual winding turns by the number of parallel paths. This is possible because each parallel circuit links the same magnetic Gux due to its spatial position in the machine. However, it should be pointed out that the winding resistance Rk (the resistance within the nite-element model) must be calculated using the parallel connection of the winding conductors. As a result, the following relationship is readily set up from Fig. 1 and Eq. (17):

dif 9A T + 2LeH {Nf } − Vf = 0; R f if + L f (18) dt 9t where Rf is the sum of all the resistive components in the eld circuit. After time discretization (using an implicit backward diHerence method), Eq. (18) can be written as:

9A 2LeH Jt T ; (19) if (t + Jt) = If (t) − {Nf } Rf Jt + Lf 9t where Jt If (t) = Rf Jt + Lf




 Lf if (t) : Vf + Jt

(20)

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

493

Fig. 2. Circuit arrangement of stator windings, generator transformer, transmission line and in nite busbar.

The nodal distribution of the vector {Nf } as follows:

eld current given by Eq. (19) is achieved using the

{Nf }if (t + Jt) = {Nf }If (t) − [Jf ]



9A 9t

eld winding

;

(21)

where [Jf ] is a densely packed square matrix (the eld winding conductivity matrix) given by [Jf ] = {Nf }{Nf }T :

(22)

Substitution of Eq. (19) into the system (14) overcomes the problem of the unknown eld current. Moreover, the eld voltage source Vf is now available to couple it with power system stabilizer and automatic voltage regulator controllers. The diHerential equations of the controllers are discretized using a similar approach and their form depends on the type of controller employed. 3.3.2. Coupling with the power network A similar procedure can be applied to the stator windings in order to connect them with the power network, represented in this work with a three-phase in nite busbar. The machine is connected to the busbar through a delta-star transformer and a transmission line represented by a series impedance. Fig. 2 shows the circuit connection arrangement. RL and LL represent the Thevenin-equivalent of the transmission line and RT and LT give the impedance per phase of the generator transformer. The stator end-winding leakage inductance is given by La , and Ra is the stator winding resistance. The objective is again to nd recurrence relations for stator currents. The following relationship can be established from the circuit of Fig. 2: [LA ]

d {I$ } + [RA ]{I$ } = −[FM ]{VI } − [M1 ]T {e}: dt

(23)

See Appendix A for explicit vector and matrix de nitions. Eq. (23) is written using a “fault” matrix [28], which allows imposing a fault condition at the high voltage terminals of the generator

494

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

transformer. Eq. (23) is similar to (18) and can be easily rearranged as  −1

1 1 T {I$ (t + Jt)} = [LA ] + [RA ] [LA ]{I$ (t)} − [FM ]{VI } − [M1 ] {e} Jt Jt

(24)

by applying the same time discretization method used for the eld circuit. The currents calculated in Eq. (24) are those circulating in the delta winding of the generator transformer, but they can be easily converted into winding currents using the [M1 ] matrix:

9A T ; (25) {IY (t + Jt)} = [M1 ]{I$ (t + Jt)} = {Is (t)} − 2LeH [TM ][N ] 9t where

−1



1 {Is (t)} = [M1 ] [LA ] + [RA ] Jt and



−1

1 [TM ] = [M1 ] [LA ] + [RA ] Jt



1 [LA ]{I$ (t)} − [FM ]{VI } Jt

[M1 ]:

(26)

(27)

Stator nodal currents are obtained by distributing the stator currents into the nite-element mesh. Thus, with the aid of the stator winding matrix [Ns ], the vector of nodal currents is given by

9A ; (28) [Ns ]{IY (t + Jt)} = [Ns ]{Is (t)} − [Js ] 9t where [Js ] is the stator conductivity matrix and is found from: [Js ] = 2LeH [Ns ][TM ][Ns ]T :

(29)

In this case, matrix [Ns ] only contains the stator winding vectors: [Ns ] = [{Na } {Nb } {Nc }]:

(30)

Once the nodal currents have been established in terms of known quantities, the space discretization of the non-linear diHusion equation (14) can be rewritten as





9A 9A 9A − [Js ] − [J] : (31) [S]{A} = {Nf }If + [Ns ]{Is } − [Jf ] 9t 9t 9t In this form, it is possible to nd a system of rst order diHerential equations with known input currents. Time discretization of the magnetic vector potential provides a system of nonlinear algebraic equations that can be solved using the well-known Newton–Raphson method [10]. Incorporation of diHerent circuit arrangements can be similarly obtained using the method outlined here. 3.4. Simulation of rotor motion The simulation of synchronous machines by nite elements requires modelling of rotor motion and several methods have been proposed to solve this problem (see for example [29–33]). An interesting

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

495

method consists in dividing the machine nite-element mesh into two parts and de ning a slip boundary at the air gap [32]. This boundary is divided into equal intervals such that nodes in the rotor side always coincide peripherally with nodes on the stator side. This means that rotor nodes belonging to the slip boundary are not allowed to fall between sides of elements in the stator side. A similar method can be devised by de ning a band of elements contained within the airgap; the element band topology being changed as the rotor moves [34], leading eventually to remeshing of the zone. Another approach is formulated using nite elements in all the machine materials but the airgap, where an analytical solution to Laplace’s equation can be found and combined with the nite element part that models the rest of the machine [33]. This permits a smooth motion of the rotor, and torque calculation is less prone to error due to mesh discretization. However, this method results in non-symmetric matrices, and increased matrix bandwidths are also obtained. Moreover, an even number of pole pitches must be used in the nite element formulation and, as a result, computer times are increased. The slip boundary method can be elegantly modi ed to enable rotor boundary nodes to fall between sides of stator elements touching the sliding interface, the magnetic vector potential of both meshes are then coupled using Lagrange multipliers [31]. This technique oHers the advantage of avoiding remeshing. However, it has been found that more iterations are needed in the solution of the system matrix when indirect methods, such as the conjugate gradient method, are utilized. Recently, Lai [29,30] has proposed a very interesting method for motion simulation where two topologically unconnected 2D nite element meshes are allowed to slide over and overlap each other, while coupled together electromagnetically. However, a complicated computational scheme must be designed to ensure that only elements of one mesh are used to model the overlapping area of the meshes. A very attractive method for modelling of rotor motion in turbine generators is devised by moving the phase bands of the stator windings with respect to the rotor and stator [23]. The stator core back is assumed to remain stationary and, hence, the small change of the reluctance path that the rotor sees in the actual machine, when movement takes place, is neglected because of the xed position of the magnetic circuit. This approach is allowable since airgaps of turbine generators are large enough to neglect tooth harmonics. The phase-band approach results in a static mesh but the time steps have to be very small to capture the system frequency oscillations. Stator winding vectors must be recalculated at each time step but the computer time required by this calculation is only a tiny fraction of the time expended in the solution of the matrix system. Furthermore, a current sheet model, where the stator phase bands are backed by a slotless stator of in nite permeability, is developed in this work to reduce computational costs; the phase band approach is naturally suited to this case because of the constant reluctance path seen by the rotor. As a result, a very eLcient and economic method for motion simulation is obtained and it is, therefore, adopted here. 3.4.1. Matrix system with motion included Once the method for motion modelling has been chosen, it is possible to solve Eq. (31) in terms of the magnetic vector potential. As it has already been explained, the phase band model requires continuous updating of the stator winding vectors and this must be reGected in the matrix system. It is now convenient to proceed with the time discretization of the magnetic vector potential, where the same time-marching method applied to winding currents is employed here. The implicit

496

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

backward diHerence procedure is unconditionally stable but small time steps are required to keep high accuracy. This is not a serious problem since the 50 Hz oscillations of stator voltages and currents need to be adequately represented, calling also for small time steps, and, therefore, errors produced by this single step-by-step method are negligible. Thus, Eq. (31) becomes:  [Jf (t + Jt)] [Js (t + Jt)] [J(t + Jt)] + + {A(t + Jt)} = {ISR }; [S(t + Jt)] + Jt Jt Jt (32) where

[Jf (t)] [Js (t)] [J(t)] + + {A(t)} {ISR } = {Nf }If (t) + [Ns (t)]{Is (t)} + Jt Jt Jt 

(33)

{ISR } can be thought of as a forcing vector, which is fully determined since it is obtained from known variables. Eq. (32) can be solved for the magnetic vector potential by inverting the matrix on the right-hand side or by applying the Newton–Raphson method if the nonlinear BH curves of the ferromagnetic materials are included. Once the solution for the time step under consideration is obtained, the next time step can be worked out in the same way as the previous one. The conductivity matrices [Jf ] and [Js ] partially disturb the sparse form of [S], as the product of two winding vectors has the eHect of connecting all the nodes of a speci c machine winding with each other. In other words, non-zero terms are placed in locations of the matrix [S] which would otherwise be empty. This results in larger matrix bandwidths and longer computing times [23]. An eLcient method can be obtained by considering products of winding vectors multiplied by the vector of magnetic vector potentials ({N }T {A}) as additional unknown variables. These variables are the Gux linkages per unit length per pole of the machine windings. Thus, Eq. (32) can be rewritten as       2LeH 2LeH [S] + [J] {N } [N ][T ] f s M {A}    {ISR }  Jt Rf Jt+Lf Jt            2LeH  T 2LeH T {Nf } {A} 0 = ; − Rf Jt+Lf {0}  Rf Jt+Lf {Nf }                2LeH 2LeH {0} t [Ns ]T {A} t+Jt [TM ]T [Ns ]T {0} [TM ]T t+Jt Jt Jt (34) where {0} is a 3 × 1 null vector. Thus, a matrix system with sparse form has been obtained at the expense of increasing the number of equations (one per machine winding) but the number of new variables is overwhelmed by the amount of nodes present in the nite-element model. Eq. (34) can then be solved using a method that exploits this sparse form. In addition, it can be seen that the matrix is always symmetric, if [TM ], which depends on the “fault matrix” [FM ], is also symmetric. This condition is commonly achieved for most of the diHerent short-circuit conditions imposed at the high voltage terminals of the generator transformer, and, therefore, the solver can also take advantage of this property. A very eLcient Gaussian elimination method is used here to solve the coupled problem. Further computer eLciency is reached by proper ordering of mesh nodes in order to bring the non-zero terms towards the main diagonal. Prediction of the magnetic vector potential at the rst iteration of

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

497

the next time step can be calculated using Euler’s method and then the calculation proceed with Eq. (34). However, the prediction stage is not strictly necessary. 3.5. The current-sheet model The computer requirements to simulate transient magnetic elds of synchronous generators by nite-element techniques are very demanding. However, it is possible to make some simpli cations to reduce the computer times demanded by the nite-element model used here for the analysis of turbine-generators operating under transient conditions. A current sheet model is proposed for this purpose, where the stator iron is assumed to have in nite permeability and the three-phase armature winding is “smeared out” and rotated with respect to the rotor, leading to a laminar distribution of stator currents. Simulation results (see below) show that the stator core back has a small inGuence, at least in machines with large air gaps. The most important result of using current sheets is that computer times are drastically reduced and that transient behavior of the solid iron rotor becomes the central area of analysis. The following considerations are made to avoid loss of accuracy: 1. The air gap length is increased to take into account the magnetomotive force drop in the stator core back. A new method, which is simple and easy to implement, was developed as follows. The eld current required to produce stator rated voltage on open circuit is calculated for the full model, then the air gap length of the current sheet model is increased until its eld current (also demanded to produce rated voltage) matches the full model current. The magnetomotive force drop in the stator is then accounted for accurately. 2. The slot leakage reactance must be added to the current-sheet model, as it is not explicitly modelled. It was calculated using simple magnetic circuit theory. The numerical value thus calculated is added to La in Fig. 2. 3. The pitch factor must be calculated because “laminar” phase bands allow for space distribution, but not for double layer short pitched windings. 4. The distribution of fractional number of winding turns to nodes using portions of element areas (see Eq. (13)) is no longer valid, as current sheets have no area. This problem is easily overcome by attributing lengths to nodes instead of areas, since the current sheets are located within an arc of element nodes. Therefore, the number of turns allocated to a particular node i can be calculated as: (i Ni = Nk ; (35) (k where (k is the angle occupied by the winding under consideration and (i is the angle attributed to node i. Once all these considerations have been made, the current sheet model is implemented in the same way as the full model. 3.6. Torque calculation The rotor speed of synchronous machines is usually not constant during transient conditions so that it is necessary to link up the mechanical system with the electromagnetic behavior of the machine.

498

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

This is done by using the swing equation which describes the rotor motion [5]: J !˙ = Tm − Te ;

(36)

where !˙ denotes the time-derivative d!=dt. Tm is the mechanical driving torque produced by the turbine. The electromagnetic torque Te can be calculated with several methods from the nite-element solution. Two important factors should be kept in mind in the selection of the torque (force) calculation method: (i) accuracy and (ii) computational cost. It has been mentioned before that transient simulations require long computer times. Therefore, the chosen method must calculate forces (torques) as quickly as possible but they must provide high accuracy at the same time. Surface integral methods meet both requirements to perform these two tasks since they only need a small part of the nite-element solution. Maxwell stresses are commonly used method of calculation in transient nite element simulations because computational costs are very low. Theoretically, Maxwell stresses must give the same value of torque (force) for any integration path enclosing the rigid body under consideration. However, nite elements give satisfactory solutions to Maxwell’s equations only in a global sense and, hence, local errors are likely to be found in the chosen integration path due to its closeness to iron-air interfaces (where local errors are mainly produced by the huge diHerences between iron and air permeabilities). On the other hand, the current sheet model lends itself to calculate the electromagnetic torque by surface integration of the Lorentz force [24]. In this case, the integration path is easily de ned because the electromagnetic torque has to be calculated from armature currents. In addition, computational costs are kept low as well. Thus, torque, in nite element terms, is given by:  Te = Bri Ii LeH rext ; (37) where Ii is the current attributed to node i and Bri is the magnetic Gux density magnitude at the ith node. Eq. (37) is seen to be the well-known BIL method [24]. The sum is carried out over all the nodes representing stator conductors and it is assumed that nodes belonging to stator windings are contained within an arc of radius rext . There is a small problem involved in Eq. (37), as Bri is not unique with rst-order elements. This is due to the magnetic Gux density discontinuities between elements that share the same node. To circumvent this problem, Bri is de ned as ni Brk ; (38) Bri = k=1 ni where Brk is the magnetic Gux density in the kth element connected to node i and the sum is performed over all the elements which have node i as a vertex. Radial magnetic Gux densities are directly worked out by expressing element shape functions in cylindrical coordinates, which provides an eLcient method of calculation since it avoids calculating the Cartesian Gux densities and from these the radial components. The driving Tm is now available to couple it with the mechanical controller in the same way as was done for the excitation controllers. 3.7. Initial condition calculation Another important aspect that must be considered in transient simulations of turbine-generators is the operating condition under which the machine is working before the transient condition. Several methods have been proposed to predict steady state operating conditions of synchronous machines under load (see for example [14,15,22]). Each of these methods is formulated under diHerent

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

499

assumptions and all of them lead to accurate prediction of eld current and load angle, but a method compatible with the phase band model is required. Refs. [14], and [15] are particularly relevant, because they use the actual distribution of stator windings and because only magnetostatic solutions are required. However, the method given in [15] demands a great deal of computations since the calculation must be done for a complete AC cycle. In addition, substantial post-processing of information is demanded. Another possibility for setting the initial condition can be found by using time-stepping techniques but they also require several solutions through time and ltering of data. Indeed, this would be the ideal procedure to obtain the initial condition, if computer times were small. Thus, it was decided that the method given in [14] is the most suitable for the initial condition calculation, as it only needs one non-linear magnetostatic solution and it is not time consuming. This technique is based on the space phasor concept which removes the complications generated by the presence of the prominent Gux third harmonic of the magnetic Gux. Appendix B shows the necessary steps to adapt this technique to the nite-element methods used in this work. 4. Results: !nite-element model validation The synchronous machine considered in this work is a two-pole 150 MVA, 120 MW, 13:8 kV, 50 Hz turbine generator. A nite-element Fortran program, that incorporates the theory given in Section 3, was developed to simulate the dynamic performance of this machine. A nite-element mesh, shown in Fig. 3, was constructed, trying to keep the number of nodes and elements as low as possible so as to reduce computational times. Nevertheless, enough elements are used in the outer surface of the rotor body to properly model the rotor eddy currents. It has 1465 nodes and 2768 elements. The transient condition simulated here is started by applying a three-phase fault at the high voltage terminals of the delta-star transformer of the synchronous machine. The generator is initially working in steady state underexcited at rated voltage, with its output power adjusted to 92% of rated output power and at 0.952 power factor. The fault is cleared and the machine is then immediately reconnected to the system after seven cycles of fault duration (0:14 ms). Test results for this machine were available and used to validate the numerical model.

Fig. 3. Finite-element mesh of the turbine generator.

500

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

Fig. 4. Steady-state magnetic Gux distribution (initial condition).

Fig. 4 shows the steady state magnetic Gux distribution calculated using the nite-element and space phasor approaches outlined in Section 3.7. The steady-state calculated load angle and eld current agree very well with test data; the errors being less than 1%. The actual machine was equipped with control systems: turbine governor and excitation system controls, whose models and parameters were also available. The excitation system uses a DC generator (controlled by an automatic voltage regulator) as a source of voltage, which provides current to the eld winding through slip rings. The governing system of the machine turbine controls the output power in the steady state and the speed under transient conditions. Figs. 5 and 6 show the comparison between test results and nite-element simulations. The nite-element model developed in this work is very eLcient, but it is important to mention that the computational times required by this model are several orders of magnitude bigger than those required by two-axis models. A time step of 0:2 ms is used to capture the 50 Hz oscillations of the machine voltage and currents. It is important to mention that the 50 Hz oscillations of the eld current, line current and output power of the test results are not shown in the gures due to a ltering process performed during the measurements. It can be seen that there is generally good agreement between test and simulated results. The 1:25 Hz oscillations of the terminal quantities following fault clearance are well predicted, although the nite element results show slightly poorer damping. The eld current, output power and load angle were somewhat underpredicted over the short circuit period, whilst the line current was overpredicted over the same period of time. The line voltage was specially well predicted in the fault period but a bit overpredicted in the instants following fault clearance. The overprediction of line current may explain the noticeable underprediction of the load angle in the rst swing, since more output power is being generated. As a result, rotor acceleration is not as high as it should be. However, subsequent oscillations show better agreement in all the machine quantities. Simulation results show a very small back swing in the load angle, which is produced by rotor losses and the interaction of the magnetic eld attached to the rotor with the eld “frozen” in stator windings at the start of the short circuit. There is also a noticeable plateau in the line voltage after fault removal, this is due to the large rotor angle which is close to 90 degrees.

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

501

Fig. 5. Comparison between test results and nite-element calculations: eld current, line current, line voltage at machine terminals.

502

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

Fig. 6. Comparison between test results and power at machine terminals.

nite-element calculations: speed deviation, rotor angle, electrical output

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

503

5. Conclusions This work presented the full development of a nite-element model of a turbine generator connected to in nite busbar through a transmission line and a delta-star connected transformer. A three-phase current sheet was proposed to avoid direct modelling of the stator coreback, leading to a reduced nite-element mesh. In addition, rotor motion is easily incorporated by simple redistribution of stator currents in the nite-element nodes that belong to the arc where the stator current sheets are located. As a result, no remeshing is necessary, avoiding the numerical problems involved with deforming meshes and the development of especial software codes demanded by sliding-surface and remeshing techniques. The loss of accuracy related to the disregarding of the stator coreback is compensated by increasing the airgap length, using a simple approach based on the open-circuit rated voltage of the machine. The electromagnetic torque is eLciently calculated by line integration of Lorentz force over the arc occupied by the three-phase current sheet. This gives a convenient method for the determination of the electromagnetic torque value, since the calculation is rapidly performed and easily programmed. Besides, the torque values is not path dependent as in the case of the Maxwell stress approach. The external devices connected to the machine windings are modelled with equivalent circuits, whose equations are strongly coupled to the nite element equations. Hence, only one non-linear solution is required per time step. The machine controls can be easily incorporated using a weak coupling. This approach has the advantage of readily changing control models without radical modi cation of the software code, which just requires the substitution of small subroutines. The nite-element techniques employed in this work were validated with test results, giving con dence in the model developed here. It is important to emphasize that the resultant nite-element model is eLcient, as it demands low computational resources. Moreover, it can be readily implemented with any programming language.

Acknowledgements The authors are grateful to the Mexican Council of Science and Technology (CONACYT) and to the National System of Researchers (SNI) of Mexico for nancial support.

Appendix A A.1. Circuit equations for the turbine-generator in;nite-busbar system The following relationships can be established from Fig. 2 by circuit analysis: {iY } = [M1 ]{i$ } = [M1 ]{IY }; {e$ } = {V } = [M1 ]{et };

504

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

{et } = −Ra [I ]{iY } − La [I ]

d{iY } − {e}; dt

{V } − {VI } = (RT + RL )[I ]{IY } + (LT + LL )[I ]

d{IY } ; dt

(A.1)

where {iY } = {ia ib ic }T

{i$ } = {iA iB iC }T ;

{e$ } = {eab ebc eca }T

{VI } = {Va Vb Vc }T ;

{et } = {eat ebt ect }T

{e} = {ea eb ec }T ;

{VI } = {VIA VIB VIC }T

{IY } = {Ia Ib Ic }T ;

(A.2)

and 

1

0

 [M1 ] =   −1

1

0

−1

−1



 0   1



1

 [I ] =  0 0

0

0



1

 0 :

0

1

(A.3)

Fault conditions at the high voltage terminals of the generator transformer are modelled by including a “fault matrix” [FM ] [28] in Eq. (A.1): {V } − [FM ]{VI } = RT [I ]{IY } + LT [I ]

d{IY } d{IY } + RL [FM ]{IY } + LL [FM ] : dt dt

The matrix [FM ] has diHerent entries according to the fault condition as follows [28]: [FM ] = [1] no fault condition.   0 0 0    [FM ] =   0 1 0  phase-a-to-earth short circuit. 0 0 1 

0

 [FM ] =  0 

0

0



0

 0  phases a and b to earth short circuit.

0

0

1

0

0

0

 [FM ] =  0 0



0

 0  three-phase-to-earth short circuit.

0

0

(A.4)

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

1  [FM ] =  

2 1 2

1 2 1 2

0

0

0

1

505



 0  phase-a-to-phase-b short circuit.

As a result, the combination of the expressions of (A.1) yields: d {I$ } + [RA ]{I$ } = −[FM ]{VI } − [M1 ]T {e}: dt where the inductance and resistance matrices [LA ] and [RA ] are given by: [LA ]

(A.5)

[LA ] = LT [I ] + LL [FM ] + La [M1 ]T [M1 ] and [RA ] = RT [I ] + RL [FM ] + Ra [M1 ]T [M1 ]: Appendix B B.1. Initial condition calculation The method of Ashtiani [14] can be set up and adapted in terms of winding vectors as follows. Fig. 7 shows the schematic diagram of one pole pitch of a two-pole synchronous machine, as well as the positions of the laminar phase bands at a speci c instant of time. The space phasors, that are required to relate the terminal conditions of the machine with its internal magnetic state, are also shown. -Ff is the magnetic Gux produced by the eld winding while -FR is the resultant magnetic Gux produced by the eld and stator windings. VF and IF are the stator terminal voltage and current.  and . are the machine load angle and power factor. The following relationships can be readily set up from Fig. 7: -Ff = -f m ej0 ; -FR = -Rm e−j ; VF = Vm e−j(+/=2) ; IF = Im e−j(+.+/=2) ;

(B.1)

where the terminal voltage VF is related to the resultant Gux -FR by: VF = −j!s -FR :

(B.2)

As a result, the resultant Gux can be written in terms of the maximum value of the terminal voltage VF and synchronous speed !s : Vm −j e : -FR = !s

(B.3)

506

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

Fig. 7. Schematic representation of space vectors of a synchronous machine.

The projection of the space phasor IF onto each phase axis gives the instantaneous values of the phase currents: ia = −im sin( + .); 
4/ ; ib = −im sin  + . − 3 
2/ ic = −im sin  + . + : 3

(B.4)

The stator Gux linkages of each phase are obtained from the nite-element magnetic vector potential solution and stator winding vectors as follows: -ih = PLeH {Ni }T {A};

i = a; b; c:

(B.5)

These Guxes are not harmonic free and they cannot be directly related to the instantaneous Guxes derived from the space projection of (B.3) onto each stator phase axis. Nevertheless, a space phasor can be constructed from the instantaneous magnetic Gux values of (B.5): 2 -F = (-ah + -bh ej2/=3 + -ch e−j2/=3 ): (B.6) 3

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

507

This equation removes the undesirable third harmonic component from the actual stator Gux linkages. The third harmonic-free Guxes are calculated by projecting this space phasor onto the stator winding axes, which is given for phase a by: -a = PLeH [ 23 {Na }T − 13 {Nb }T − 13 {Nc }T ]{A} = PLeH {Nafil }T {A};

(B.7)

where {Nafil } is a ‘new’ winding vector that has the additional property of eliminating the presence of the third harmonic Gux component. Similar winding vectors can be readily established for phases b and c. These harmonic free Guxes can be now related to the Gux projection of (B.3) as follows: Vm cos  = PLeH {Nafil }T {A}; !s 
Vm 2/ = PLeH {Nbfil }T {A}: cos  + -b = !s 3 -a =

(B.8)

Only two Gux equations are required since the Gux linkages of phase c can be derived from (B.8). Eqs. (B.4) and (B.8) can be combined with the steady-state nite-element model of the synchronous machine (Eq. (31) without the term [J]{9A=9t}) in order to obtain the following set of equations: {E1 } = [S]{A} − {Nf }if − {Na }ia − {Nb }ib − {Nc }ic ; Vm cos ; !s 
Vm 2/ T ; cos  + E3 = PLeH {Nbfil } {A} − !s 3 E2 = PLeH {Nafil }T {A} −

(B.9)

where {E1 } is a null vector when the solution is found. Similarly, E2 and E1 are required to be zero at the solution point. This system of non-linear equations is iteratively solved using the Newton– Rhapson method [10]. This procedure produces a $ increment vector which is used to obtain an improved solution at the (k + 1)th iteration. This increment vector is given by: k+1 k   k  [H ] {I  } −{Nf }   {JA}  {E1 }            T V m  PLeH {Nafil }  sin  0 J E = − (B.10) 2 !s            Ji   E  ) 0 PLeH {Nbfil }T V!ms sin( + 2/ f 3 3 where [H ] is the Hessian matrix of [S]{A} and it is calculated as: [H ] =

9 [S]{A} 9{A}

(B.11)

{I  } is a vector of currents that is obtained by diHerentiating (B.4) and distributing the resulting currents into their respective phases:   

4/ 2/  {I } = Im {Na }cos( + .) + {Nb } cos  + . − + {Nc } cos  + . − : (B.12) 3 3

508

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

References [1] R.H. Park, Two-reaction theory of synchronous machines. Generalized method of analysis, AIEE Trans. 48 (1929) 716–730. [2] I.M. Canay, Causes of discrepancies on calculation of rotor quantities and exact equivalent diagrams of the synchronous machine, IEEE Trans. Power Apparatus Syst. 88 (1969) 1114–1120. [3] I.M. Canay, Modelling of alternating-current machines having multiple rotor circuits, IEEE Trans. Energy Conversion 8 (1993) 280–296. [4] IEEE Guide: Test Procedures for synchronous machines, IEEE Std. 115-1995. [5] P. Kundur, Power System Stability and Control, McGraw Hill, New York, 1994. [6] I. Kamwa, P. Viarouge, H. Le-Huy, J. Dickinson, A frequency domain maximum likelihood estimation of synchronous machine high-order models using SSFR test data, IEEE Trans. Energy Conversion 7 (1992) 525–536. [7] A. Keyhani, H. Tsai, Identi cation of high order synchronous generator models from SSFR test data, IEEE Trans. Energy Conversion 9 (1994) 593–603. [8] R. Escarela-Perez, T. Niewierowicz, E. Campero-Littlewood, Synchronous machine parameters from frequencyresponse nite-element simulations and genetic algorithms, IEEE Trans. Energy Conversion 16 (2001) 198–203. [9] T. Niewierowicz, R. Escarela-Perez, E. Campero-Littlewood, Hybrid genetic algorithm for the identi cation of high-order synchronous machine two-axis equivalent circuits, Inter. J. Comput. Electrical Eng. 29 (4) (2003) 505–522. [10] P.P. Silvester, R.L. Ferrari, Finite Elements for Electrical Engineers, 3rd Edition, Cambridge University Press, Cambridge, 1996. [11] M.V.K. Chari, P. Silvester, Analysis of turboalternator magnetic elds by nite elements, IEEE Trans. Power Apparatus Syst. 90 (1971) 454–464. [12] P. Brandl, K. Reichert, W. Vogt, Simulation of turbogenerators on steady state load, Brown Bovery Rev. 9 (1975) 444 – 449. [13] M.V.K. Chari, S.H. Minnich, S.C. Tandon, Z.J. Csendes, J. Berkery, Load characteristics of synchronous generators by the nite element method, IEEE Trans. Power Apparatus Syst. 100 (1981) 1–9. [14] C.N. Ashtiani, D.A. Lowther, The use of nite elements in the simulation of steady state operation of a synchronous generator with a known terminal loading condition, IEEE Trans. Magnetics 19 (1983) 2381–2384. [15] S.R. Chaudhry, S. Ahmed-Zaid, N.A. Demerdash, Coupled nite-element/state-space modeling of turbogenerators in the ABC frame of reference—the short-circuit and load cases including saturated parameters, IEEE Trans. Energy Conversion 10 (1995) 63–70. [16] A.Y. Hannalla, D.C. Macdonald, Numerical analysis of transient eld problems in electrical machines, IEE Proc. (Pt. C) 123 (1976) 893–898. [17] A.Y. Hannalla, D.C. Macdonald, Sudden 3-phase short-circuit characteristics of turbine-generator from design data using electromagnetic eld calculations, IEE Proc. (Pt. C) 127 (1980) 213–220. [18] J.P. Sturgess, M. Zhu, D.C. Macdonald, Finite-element simulation of a generator on load during and after a three-phase fault, IEEE Trans. Energy Conversion 7 (1992) 787–793. [19] S.J. Salon, M.J. DeBortoli, R. Palma, Coupling of transient elds, circuits and motion using nite element analysis, J. Electromagnetic Waves Appl. 4 (1990) 1077–1106. [20] G. Bedrosian, A new method for coupling nite element eld solution with external circuits and kinematics, IEEE Trans. Magnetics 29 (1993) 1664–1668. [21] I. Tsukerman, Accurate computation of ‘ripple solutions’ on moving nite element meshes, IEEE Trans. Magnetics 31 (1995) 1472–1475. [22] P.J. Turner, D.C. Macdonald, Finite element prediction of steady state operation of a turbine-generator using conductivity matrices, IEEE Trans. Magnetics 17 (1981) 3262–3264. [23] P.J. Turner, Finite-element simulation of turbine-generator terminal faults and application to machine parameter prediction, IEEE Trans. Energy Conversion 2 (1987) 122–131. [24] W.H. Hayt Jr., Engineering Electromagnetics, 5th edition, McGraw Hill, New York, 1989. [25] R. Escarela-Perez, D.C. Macdonald, A novel nite-element transient computation of two-axis parameters of solid-rotor generators for use in power systems, IEEE Trans. Energy Conversion 13 (1998) 49–54.

R. Escarela-Perez et al. / Finite Elements in Analysis and Design 40 (2004) 485 – 509

509

[26] J.R. Brauer, B.E. MacNeal, L.A. Larkin, V.D. Overbye, New method of modelling electronic circuits coupled with 3D electromagnetic nite element models, IEEE Trans. Magnetics 27 (1991) 4085–4088. [27] S.L. Ho, L. Li, W.N. Fu, H.C. Wong, A novel approach to circuit- eld-torque coupled time stepping nite element modeling of electric machines, IEEE Trans. Magentics 36 (2000) 1886–1889. [28] M.A. Barber, M. Giannini, Simulation of a synchronous generator connected via a delta-star transformer, IEE Proc. 121 (1974) 1513–1521. [29] H.C. Lai, P.C. Coles, D. Rodger, P.J. Leonard, Transient analysis of an electromagnetic actuator using an overlapping nite element scheme, IEEE Trans. Magnetics 36 (2000) 1462–1467. [30] H.C. Lai, D. Rodger, P.C. Coles, A nite element scheme for colliding meshes, IEEE Trans. Magnetics 35 (1999) 1362–1364. [31] D. Rodger, H.C. Lai, P.J. Leonard, Coupled elements for problems involving movement, IEEE Trans. Magnetics 26 (1990) 548–550. [32] T.W. Preston, A.B.J. Reece, P.S. Sangha, Induction motor analysis by time-stepping techniques, IEEE Trans. Magnetics 24 (1988) 471–474. [33] A.A. Abdel-Razek, J.L. Coulomb, M. Feliachi, J.C. Sabonnadiere, Conception of an air-gap element for the dynamic analysis of the electromagnetic eld in electric machines, IEEE Trans. Magnetics 18 (1982) 655–659. [34] B. Davat, Z. Ren, Lajoie-Mazenc, The movement in eld modeling, IEEE Trans. Magnetics 21 (1985) 2296 –2298.