Massachusetts Institute of Technology

Department of Electrical Engineering and Computer Science

6.685 Electric Machinery

Class Notes 9: Synchronous Machine Simulation Models September 5, 2005

c
�2005 James L. Kirtley Jr.

1 Introduction
In this document we develop models useful for calculating the dynamic behavior of synchronous
machines. We start with a commonly accepted picture of the synchronous machine, assuming that
the rotor can be fairly represented by three equivalent windings: one being the field and the other
two, the d- and q- axis “damper” windings, representing the effects of rotor body, wedge chain,
amortisseur and other current carrying paths.
While a synchronous machine is assumed here, the results are fairly directly applicable to
induction machines. Also, extension to situations in which the rotor representation must have
more than one extra equivalent winding per axis should be straightforward.

2 Phase Variable Model

To begin, assume that the synchronous machine can be properly represented by six equivalent
windings. Four of these, the three armature phase windings and the field winding, really are
windings. The other two, representing the effects of distributed currents on the rotor, are referred
to as the “damper” windings. Fluxes are, in terms of currents:
" # " #" #
λph Lph M Iph
= (1)
λR MT LR IR

where phase and rotor fluxes (and, similarly, currents) are:

 
λa
λph =  λb  (2)
 
λc
 
λf
λR =  λkd  (3)
 
λkq
There are three inductance sub- matrices. The first of these describes armature winding induc-
tances:  
La Lab Lac
Lph =  Lab Lb Lbc  (4)
 
Lac Lbc Lc

where, for a machine that may have some saliency:

La = La0 + L2 cos 2θ (5)
2π
Lb = La0 + L2 cos 2(θ − ) (6)
3
2π
Lc = La0 + L2 cos 2(θ + ) (7)
3
π
Lab = Lab0 + L2 cos 2(θ − ) (8)
3
Lbc = Lab0 + L2 cos 2θ (9)
π
Lac = Lab0 + L2 cos 2(θ + ) (10)
3
Note that, in this last set of expressions, we have assumed a particular form for the mutual in-
ductances. This is seemingly restrictive, because it constrains the form of phase- to- phase mutual
inductance variations with rotor position. The coefficient L2 is actually the same in all six of these
last expressions. As it turns out, this assumption does not really restrict the accuracy of the model
very much. We will have more to say about this a bit later.
The rotor inductances are relatively simply stated:
 
Lf Lf kd 0
LR =  Lf kd Lkd 0  (11)
 
0 0 Lkq
And the stator- to- rotor mutual inductances are:
 
M cos θ Lakd cos θ −Lakq sin θ
2π 2π 2π
M =  M cos(θ − 3 ) Lakd cos(θ − 3 ) −Lakq sin(θ − 3 ) (12)
 

2π 2π 2π
M cos(θ + 3 ) Lakd cos(θ + 3 ) −Lakq sin(θ + 3 )

3 Park’s Equations
The first step in the development of a suitable model is to transform the armature winding variables
to a coordinate system in which the rotor is stationary. We identify equivalent armature windings
in the direct and quadrature axes. The direct axis armature winding is the equivalent of one of
the phase windings, but aligned directly with the field. The quadrature winding is situated so
that its axis leads the field winding by 90 electrical degrees. The transformation used to map
the armature currents, fluxes and so forth onto the direct and quadrature axes is the celebrated
Park’s Transformation, named after Robert H. Park, an early investigator into transient behavior
in synchronous machines. The mapping takes the form:
   
ud ua
u = u = T u = T  ub  (13)
   
 q  dq ph
u0 uc
Where the transformation and its inverse are:
 
cos θ cos(θ − 23π ) cos(θ + 23π )
2
T =  − sin θ − sin(θ − 3 ) − sin(θ + 23π ) 
2π
(14)

3 1 1 1
2 2 2

2
  
cos θ − sin θ 1
T −1 =  cos(θ − 2π ) − sin(θ − 2π
) 1  (15)
 
3 3
2π 2π
cos(θ + 3 ) − sin(θ + 3 ) 1
This transformation maps balanced sets of phase currents into constant currents in the d-q frame.
That is, if rotor angle is θ = ωt + θ0 , and phase currents are:

Ia = I cos ωt
2π
Ib = I cos(ωt − )
3
2π
Ic = I cos(ωt + )
3
Then the transformed set of currents is:

Id = I cos θ0
Iq = −I sin θ0

Now, we apply this transformation to (1) to express fluxes and currents in the armature in the d-q
reference frame. To do this, extract the top line in (1):

λph = Lph I ph + M I R (16)

The transformed flux is obtained by premultiplying this whole expression by the transformation
matrix. Phase current may be obtained from d-q current by multiplying by the inverse of the
transformation matrix. Thus:
λdq = T Lph T −1 I dq + T M I R (17)
The same process carried out for the lower line of (1) yields:

λR = M T T −1 I dq + LR I R (18)

Thus the fully transformed version of (1) is:

" # " #" #
λdq Ldq LC Idq
= 3 T (19)
λR 2 LC LR IR

If the conditions of (5) through (10) are satisfied, the inductance submatrices of (19) wind up being
of particularly simple form. (Please note that a substantial amount of algebra has been left out
here!)  
Ld 0 0
Ldq =  0 Lq 0  (20)
 
0 0 L0
 
M Lakd 0
LC =  0 0 Lakq  (21)
 
0 0 0

3
 Note that (19) through (21) express three separate sets of apparently independent flux/current
relationships. These may be re-cast into the following form:
    
λd Ld Lakd M Id
3
 λkd  =  L L L Ikd  (22)
    
2 akd kd f kd  
3
λf 2M Lf kd Lf If
" # " #" #
λq Lq Lakq Iq
= 3 (23)
λkq L
2 akq Lkq Ikq

λ0 = L0 I0 (24)
Where the component inductances are:

3
Ld = La0 − Lab0 + L2 (25)
2
3
Lq = La0 − Lab0 − L2 (26)
2
L0 = La0 + 2Lab0 (27)

Note that the apparently restrictive assumptions embedded in (5) through (10) have resulted in
the very simple form of (21) through (24). In particular, we have three mutually independent sets
of fluxes and currents. While we may be concerned about the restrictiveness of these expressions,
note that the orthogonality between the d- and q- axes is not unreasonable. In fact, because these
axes are orthogonal in space, it seems reasonable that they should not have mutual flux linkages.
The principal consequence of these assumptions is the de-coupling of the zero-sequence component
of flux from the d- and q- axis components. We are not in a position at this time to determine
the reasonableness of this. However, it should be noted that departures from this form (that is,
coupling between the “direct” and “zero” axes) must be through higher harmonic fields that will
not couple well to the armature, so that any such coupling will be weak.
Next, armature voltage is, ignoring resistance, given by:

d d
λ = T −1 λdq
V ph = (28)
dt ph dt
and that the transformed armature voltage must be:

V dq = T V ph
d
= T (T −1 λdq )
dt
d d
= λdq + (T T −1 )λdq (29)
dt dt
A good deal of manupulation goes into reducing the second term of this, resulting in:
 
0 − dθ 0
d −1  dθ
dt
T T = 0 0  (30)

dt dt
0 0 0

4
This expresses the speed voltagethat arises from a coordinate transformation. The two voltage/flux
relationships that are affected are:
dλd
Vd = − ωλq (31)
dt
dλq
Vq = + ωλd (32)
dt
where we have used
dθ
ω= (33)
dt

4 Power and Torque

Instantaneous power is given by:
P = V a Ia + V b Ib + V c Ic (34)
Using the transformations given above, this can be shown to be:
3 3
P = Vd Id + Vq Iq + 3V0 I0 (35)
2 2
which, in turn, is:

3 3 dλd dλq dλ0

P = ω (λd Iq − λq Id ) + ( Id + Iq ) + 3 I0 (36)
2 2 dt dt dt
Then, noting that electrical speed ω and shaft speed Ω are related by ω = pΩ and that (36)
describes electrical terminal power as the sum of shaft power and rate of change of stored energy,
we may deduce that torque is given by:
3
T = p(λd Iq − λq Id ) (37)
2

5 Per-Unit Normalization
The next thing for us to do is to investigate the way in which electric machine system are nor-
malized, or put into what is called a per-unit system. The reason for this step is that, when the
voltage, current, power and impedance are referred to normal operating parameters, the behavior
characteristics of all types of machines become quite similar, giving us a better way of relating
how a particular machine works to some reasonable standard. There are also numerical reasons for
normalizing performance parameters to some standard.
The first step in normalization is to establish a set of base quantities. We will be normalizing
voltage, current, flux, power, impedance and torque, so we will need base quantities for each of
these. Note, however, that the base quantities are not independent. In fact, for the armature, we
need only specify three quantities: voltage (VB ), current (IB ) and frequency (ω0 ). Note that we do
not normalize time nor frequency. Having done this for the armature circuits, we can derive each
of the other base quantities:

5
 • Base Power
3
PB = VB IB
2
• Base Impedance
VB
ZB =
IB
• Base Flux
VB
λB =
ω0
• Base Torque
p
TB = PB
ω0

Note that, for our purposes, base voltage and current are expressed as peak quantities. Base voltage
is taken on a phase basis (line to neutral for a “wye” connected machine), and base current is
similarly taken on a phase basis, (line current for a “wye” connected machine).
Normalized, or per-unit quantities are derived by dividing the ordinary variable (with units) by
the corresponding base. For example, per-unit flux is:

λ ω0 λ
ψ= = (38)
λB VB
In this derivation, per- unit quantities will usually be designated by lower case letters. Two
notable exceptions are flux, where we use the letter ψ, and torque, where we will still use the upper
case T and risk confusion.
Now, we note that there will be base quantities for voltage, current and frequency for each of
the different coils represented in our model. While it is reasonable to expect that the frequency
base will be the same for all coils in a problem, the voltage and current bases may be different. We
might write (22) as:
 ω0 IdB ω0 IkB ω0 I f B 
Vdb Ld Vdb Lakd Vdb M
   
ψd   id
ω0 IdB 3 ω0 IkB ω0 I f B
 ψkd  =  Vkb 2 Lakd Vkb Lkd Vkdb Lf kd   ikd  (39)
    
ψf if
 
ω0 IdB 3 ω0 IkB ω0 I f B
Vf b 2 M Vf b Lf kd Vf b Lf

where i = I/IB denotes per-unit, or normalized current.

Note that (39) may be written in simple form:
    
ψd xd xakd xad id
 ψkd  =  xakd xkd xf kd   ikd  (40)
    
ψf xad xf kd xf if

It is important to note that (40) assumes reciprocity in the normalized system. To wit, the following
expressions are implied:
IdB
xd = ω 0 Ld (41)
VdB

6
 IkB
xkd = ω0 Lkd (42)
VkB
If B
xf = ω0 Lf (43)
Vf B
IkB
xakd = ω0 Lakd
VdB
3 IdB
= ω0 Lakd (44)
2 VkB
If B
xad = ω0 M
VdB
3 IdB
= ω0 M (45)
2 Vf B
IkB
xf kd = ω0 Lf kd
Vf b
If B
= ω0 Lf kd (46)
Vkb
These in turn imply:

3
VdB IdB = Vf B If B (47)
2
3
VdB IdB = VkB IkB (48)
2
Vf B If B = VkB IkB (49)
These expressions imply the same power base on all of the windings of the machine. This is
so because the armature base quantities Vdb and Idb are stated as peak values, while the rotor base
quantities are stated as DC values. Thus power base for the three- phase armature is 23 times
the product of peak quantities, while the power base for the rotor is simply the product of those
quantities.
The quadrature axis, which may have fewer equivalent elements than the direct axis and which
may have different numerical values, still yields a similar structure. Without going through the
details, we can see that the per-unit flux/current relationship for the q- axis is:
" # " #" #
ψq xq xakq iq
= (50)
ψkq xakq xkq ikq

The voltage equations, including speed voltage terms, (31) and (32), may be augmented to
reflect armature resistance:
dλd
Vd = − ωλq + Ra Id (51)
dt
dλq
Vq = ωλd + + Ra Iq (52)
dt
The per-unit equivalents of these are:
1 dψd ω
vd = − ψq + ra id (53)
ω0 dt ω0

 ω 1 dψq
vq = ψd + + ra iq (54)
ω0 ω0 dt

Where the per-unit armature resistance is just ra = ZRBa

Note that none of the other circuits in this model have speed voltage terms, so their voltage
expressions are exactly what we might expect:
1 dψf
vf = + rf if (55)
ω0 dt
1 dψkd
vkd = + rkd ikd (56)
ω0 dt
1 dψkq
vkq = + rkq ikq (57)
ω0 dt
1 dψ0
v0 = + ra i0 (58)
ω0 dt
It should be noted that the damper winding circuits represent closed conducting paths on the rotor,
so the two voltages vkd and vkq are always zero.
Per-unit torque is simply:
Te = ψd iq − ψq id (59)
Often, we need to represent the dynamic behavior of the machine, including electromechanical
dynamics involving rotor inertia. If we note J as the rotational inertia constant of the machine
system, the rotor dynamics are described by the two ordinary differential equations:
1 dω
J = T e + T m (60)
p dt
dδ
= ω − ω0 (61)
dt
where T e and T m represent electrical and mechanical torques in “ordinary” variables. The angle δ
represents rotor phase angle with respect to some synchronous reference.
It is customary to define an “inertia constant” which is not dimensionless but which nevertheless
fits into the per-unit system of analysis. This is:
Rotational kinetic energy at rated speed
H ≡ (62)
Base Power
Or:  2
1 ω0
2J p J ω0
H= = (63)
PB 2pTB
Then the per-unit equivalent to (60) is:

2H dω
= Te + Tm (64)
ω0 dt
where now we use Te and Tm to represent per-unit torques.

8
6 Equal Mutual’s Base
In normalizing the differential equations that make up our model, we have used a number of base
quantities. For example, in deriving (40), the per-unit flux- current relationship for the direct
axis, we used six base quantities: VB , IB , Vf B , If B , VkB and IkB . Imposing reciprocity on (40)
results in two constraints on these six variables, expressed in (47) through (49). Presumably the
two armature base quantities will be fixed by machine rating. That leaves two more “degrees of
freedom” in selection of base quantities. Note that the selection of base quantities will affect the
reactance matrix in (40).
While there are different schools of thought on just how to handle these degrees of freedom, a
commonly used convention is to employ what is called the equal mutuals base system. The two
degrees of freedom are used to set the field and damper base impedances so that all three mutual
inductances of (40) are equal:
xakd = xf kd = xad (65)
The direct- axis flux- current relationship becomes:
    
ψd xd xad xad id
 ψkd  =  xad xkd xad   ikd  (66)
    
ψf xad xad xf if

7 Equivalent Circuit

i�
d ra xal xf l if
rf �
∧∧∧
∨∨
∩∩∩∩ ∩∩∩∩ ∧∧∧
∨∨
+ + ⊃ +
⊃ xkdl
⊃
⊃ ⊃

(ω0 vd + ωψq ) ψd ⊃ xad vf

⊃
⊃
<>
<> rkd
- - < -

Figure 1: D- Axis Equivalent Circuit

The flux- current relationship of (66) is represented by the equivalent circuit of Figure 1, if the
“leakage” inductances are defined to be:

xal = xd − xad (67)

xkdl = xkd − xad (68)
xf l = xf − xad (69)

9
 Many of the interesting features of the electrical dynamics of the synchronous machine may be
discerned from this circuit. While a complete explication of this thing is beyond the scope of this
note, it is possible to make a few observations.
The apparent inductance measured from the terminals of this equivalent circuit (ignoring resis-
tance ra ) will, in the frequency domain, be of the form:

ψd (s) Pn (s)
x(s) = = xd (70)
id (s) Pd (s)

Both the numerator and denominator polynomials in s will be second order. (You may convince
yourself of this by writing an expression for terminal impedance). Since this is a “diffusion” type
circuit, having only resistances and inductances, all poles and zeros must be on the negative real
axis of the “s-plane”. The per-unit inductance is, then:

(1 + Td′ s)(1 + Td′′ s)

x(s) = xd ′ s)(1 + T ′′ s) (71)
(1 + Tdo do
The two time constants Td′ and Td′′ are the reciprocals of the zeros of the impedance, which are
the poles of the admittance. These are called the short circuittime constants.
The other two time constants Tdo′ and T ′′ are the reciprocals of the poles of the impedance, and
do
so are called the open circuittime constants.
We have cast this thing as if there are two sets of well- defined time constants. These are the
transient time constants Td′ and Tdo′ , and the subtransient time constants T ′′ and T ′′ . In many
d do
cases, these are indeed well separated, meaning that:

Td′ ≫ Td′′ (72)

′ ′′
Tdo ≫ Tdo (73)

If this is true, then the reactance is described by the pole-zero diagram shown in Figure 2.
Under this circumstance, the apparent terminal inductance has three distinct values, depending on
frequency. These are the synchronous inductance, the transient inductance, and the subtransient
inductance, given by:

Td′
x′d = xd ′ (74)
Tdo
T ′′
x′′d = xd′ d′′
Tdo
T ′ T ′′
= xd ′d d′′ (75)
Tdo Tdo

A Bode Plotof the terminal reactance is shown in Figure 3.

If the time constants are spread widely apart, they are given, approximately, by:

′ xf
Tdo = (76)
ω0 rf
′′ xkdl + xf l ||xad
Tdo = (77)
ω0 rkd

10
 1 1
′

Tdo
� �
Tdo ”
× ×
1 1
Td ” Td′

Figure 2: Pole-Zero Diagram For Terminal Inductance

log |x(jω)|
�
��
�
��

1
′
1 1 1 log ω
Tdo T ′
d Tdo ” Td ”

Figure 3: Frequency Response of Terminal Inductance

Finally, note that the three reactances are found simply from the model:

xd = xal + xad (78)

x′d = xal + xad ||xf l (79)
x′′d = xal + xad ||xf l ||xkdl (80)

8 Statement of Simulation Model

Now we can write down the simulation model. Actually, we will derive more than one of these, since
the machine can be driven by either voltages or currents. Further, the expressions for permanent
magnet machines are a bit different. So the first model is one in which the terminals are all
constrained by voltage.
The state variables are the two stator fluxes ψd , ψq , two “damper” fluxes ψkd , ψkq , field flux ψf ,
and rotor speed ω and torque angle δ. The most straightforward way of stating the model employs
currents as auxiliary variables, and these are:
   −1  
id xd xad xad ψd
 ikd  =  xad xkd xad   ψkd  (81)
     
if xad xad xf ψf
" # " #−1 " #
iq xq xaq ψq
= (82)
ikq xaq xkq ψkq

11
Then the state equations are:
dψd
= ω0 vd + ωψq − ω0 ra id (83)
dt
dψq
= ω0 vq − ωψd − ω0 ra iq (84)
dt
dψkd
= −ω0 rkd ikd (85)
dt
dψkq
= −ω0 rkq ikq (86)
dt
dψf
= ω0 vf − ω0 rf if (87)
dt
dω ω0
= (Te + Tm ) (88)
dt 2H
dδ
= ω − ω0 (89)
dt
and, of course,
Te = ψd iq − ψq id

8.1 Statement of Parameters:

Note that often data for a machine may be given in terms of the reactances xd , x′d , xd′′ , Tdo
′ and
′′ , rather than the elements of the equivalent circuit model. Note that there are four inductances
Tdo
in the equivalent circuit so we have to assume one. There is no loss in generality in doing so.
Usually one assumes a value for the stator leakage inductance, and if this is done the translation is
straightforward:

xad = xd − xal
xad (x′d − xal )
xf l =
xad − x′d + xal
1
xkdl = 1 1 1
x′′
−xal − xad −
d
xf l

xf l + xad
rf = ′
ω0 Tdo
xkdl + xad ||xf l
rkd = ′′
ω0 Tdo

8.2 Linearized Model

Often it becomes desirable to carry out a linearized analysis of machine operation to, for example,
examine the damping of the swing mode at a particular operating point. What is done, then,
is to assume a steady state operating point and examine the dynamics for deviations from that
operating point that are “small”. The definition of “small” is really “small enough” that everything
important appears in the first-order term of a Taylor series about the steady operating point.
Note that the expressions in the machine model are, for the most part, linear. There are,
however, a few cases in which products of state variables cause us to do the expansion of the

12

Taylor series. Assuming a steady state operating point [ψd0 ψkd0 ψf 0 ψq0 ψkq0 ω0 δ0 ], the first-
order (small-signal) variations are described by the following set of equations. First, since the
flux-current relationship is linear:
   −1  
id1 xd xad xad ψd1
 ikd1  =  xad xkd xad   ψkd1  (90)
     
if 1 xad xad xf ψf 1
" # " #−1 " #
iq1 xq xaq ψq1
=
(91)
ikq1 xaq xkq ψkq1

Terminal voltage will be, for operation against a voltage source:

Vd = V sin δ Vq = V cos δ

Then the differential equations governing the first-order variations are:

dψd1
= ω0 V cos δ0 δ1 + ω0 ψq1 + ω1 ψq0 − ω0 ra id1 (92)
dt

dψq1

= −ω0 V sin δ0 δ1 − ω0 ψd1 − ω1 ψd0 − ω0 ra iq1 (93)

dt

dψkd1

= −ω0 rkd ikd1 (94)

dt

dψkq1

= −ω0 rkq ikq1 (95)

dt

dψf 1

= −ω0 rf if 1 (96)
dt
dω1 ω0
= (Te1 + Tm1 ) (97)
dt 2H

dδ1

= ω1 (98)
dt
Te = ψd0 iq1 + ψd1 iq0 − ψq0 id1 − ψq1 id0

8.3 Reduced Order Model for Electromechanical Transients

In many situations the two armature variables contribute little to the dynamic response of the
machine. Typically the armature resistance is small enough that there is very little voltage drop
across it and transients in the difference between armature flux and the flux that would exist in
the “steady state” decay rapidly (or are not even excited). Further, the relatively short armature
time constant makes for very short time steps. For this reason it is often convenient, particularly
when studying the relatively slow electromechanical transients, to omit the first two differential
equations and set:

ψd = vq = V cos δ (99)
ψq = −vd = −V sin δ (100)

The set of differential equations changes only a little when this approximation is made. Note,
however, that it can be simulated with far fewer “cycles” if the armature time constant is short.

13

9 Current Driven Model: Connection to a System

The simulation expressions developed so far are useful in a variety of circumstances. They are,
however, difficult to tie to network simulation programs because they use terminal voltage as an
input. Generally, it is more convenient to use current as the input to the machine simulation and
accept voltage as the output. Further, it is difficult to handle unbalanced situations with this set
of equations.
An alternative to this set would be to employ the phase currents as state variables. Effectively,
this replaces ψd , ψq and ψ0 with ia , ib , and ic . The resulting model will, as we will show, interface
nicely with network simulations.
To start, note that we could write an expression for terminal flux, on the d- axis:
xad ||xkdl xad ||xf l
ψd = x′′d id + ψf + ψkd (101)
xad ||xkdl + xf l xad ||xf l + xkdl
and here, of course,
x′′d = xal + xad ||xkdl ||xf l
This leads us to define a “flux behind subtransient reactance”:
xad xkdl ψf + xad xf l ψkd
ψd′′ = (102)
xad xkdl + xad xf l + xkdl xf l
So that
ψd = ψd′′ + xd′′ id
On the quadrature axis the situation is essentially the same, but one step easier if there is only
one quadrature axis rotor winding:
xaq
ψq = x′′q iq + ψkq (103)
xaq + xkql
where
x′′q = xal + xaq ||xkql
Very often these fluxes are referred to as “voltage behind subtransient reactance, with ψd′′ = eq′′
and ψq′′ = −ed′′ . Then:

ψd = x′′d id + e′′q (104)

ψq = x′′q iq − ed′′ (105)

Now, if id and iq are determined, it is a bit easier to find the other currents required in the
simulation. Note we can write:
" # " #" # " #
ψkd xkd xad ikd xad
= + id (106)
ψf xad xf if xad

and this inverts easily:

" # " #−1 " # " # !
ikd xkd xad ψkd xad
= − id (107)
if xad xf ψf xad

14
The quadrature axis rotor current is simply:
1 xaq
ikq = ψkq − iq (108)
xkq xkq

The torque equation is the same, but since it is usually convenient to assemble the fluxes behind
subtransient reactance, it is possible to use:

Te = e′′q iq + ed′′ id + (xd′′ − xq′′ )id iq (109)

Now it is necessary to consider terminal voltage. This is most conveniently cast in matrix
notation. The vector of phase voltages is:
 
va
v ph =  vb  (110)
 
vc
Then, with similar notation for phase flux, terminal voltage is, ignoring armature resistance:

1 dψ ph
v ph =
ω0 dt
1 d n −1 o
= T ψ dq (111)
ω0 dt
Note that we may define the transformed vector of fluxes to be:

ψ dq = x′′ idq + e′′ (112)

where the matrix of reactances shows orthogonality:
 
x′′ 0 0
′′  d
x =  0 x′′q 0  (113)

0 0 x0
and the vector of internal fluxes is:
 
e′′q
′′
e =  −e ′′d  (114)
 
0
Now, of course, idq = T iph , so that we may re-cast (111) as:

1 d n −1 ′′ o
v ph = T x T iph + T −1 e ′′ (115)
ω0 dt
Now it is necessary to make one assumption and one definition. The assumption, which is
only moderately restrictive, is that subtransient saliency may be ignored. That is, we assume
that x′′d = x′′q . The definition separates the “zero sequence” impedance into phase and neutral
components:

15
 x0 = x′′d + 3xg (116)
Note that according to this definition the reactance xg accounts for any impedance in the neutral
of the synchronous machine as well as mutual coupling between phases.
Then, the impedance matrix becomes:
   
x′′ 0 0 0 0 0
′′  d
x =  0 x′′d 0  +  0 0 0  (117)
  
0 0 x′′d 0 0 3xg
In compact notation, this is:

x′′ = x′′d I + xg (118)

where I is the identity matrix.
Now the vector of phase voltages is:
1 d n ′′ o
v ph = xd iph + T −1 xg T iph + T −1 e′′ (119)
ω0 dt
Note that in (119), we have already factored out the multiplication by the identity matrix. The
next step is to carry out the matrix multiplication in the third term of (119). This operation turns
out to produce a remarkably simple result:
 
1 1 1
−1
T xg T = xg  1 1 1  (120)
 
1 1 1
The impact of this is that each of the three phase voltages has the same term, and that is
related to the time derivative of the sum of the three currents, multiplied by xg .
The third and final term in (119) describes voltages induced by rotor fluxes. It can be written
as:

1 d n −1 ′′ o 1 d n −1 o ′′ 1 de′′
T e = T e + T −1 (121)
ω0 dt ω0 dt ω0 dt
Now, the time derivative of the inverse transform is:
 
− sin(θ) − cos(θ) 0
1 d −1 ω  2π 2π
T = − sin(θ − ) − cos(θ − ) 0  (122)

ω0 dt ω0
 3 3
− sin(θ + 23π ) − cos(θ + 2π
3 ) 0
Now the three phase voltages can be extracted from all of this matrix algebra:

x′′d dia xg d
va = + (ia + ib + ic ) + e′′a (123)
ω0 dt ω0 dt
x′′d dib xg d
vb = + (ia + ib + ic ) + eb′′ (124)
ω0 dt ω0 dt
x′′d dic xg d
vc = + (ia + ib + ic ) + ec′′ (125)
ω0 dt ω0 dt

16
 Where the internal voltages are:

ω ′′
e′′a = − (e sin(θ) − ed′′ cos(θ))
ω0 q
1 de′′q 1 de′′
+ cos(θ) + sin(θ) d (126)
ω0 dt ω0 dt
ω ′′ 2π 2π
e′′b = − (eq sin(θ − ) − e′′d cos(θ − ))
ω0 3 3
1 2π de′′q 1 2π de′′d
+ cos(θ − ) + sin(θ − ) (127)
ω0 3 dt ω0 3 dt
ω 2π 2π
e′′c = − (e′′q sin(θ + ) − e′′d cos(θ + ))
ω0 3 3
1 2π de′′q 1 2π de′′d
+ cos(θ + ) + sin(θ + ) (128)
ω0 3 dt ω0 3 dt
This set of expressions describes the equivalent circuit shown in Figure 4.

′′
i� ′′
a xd ��
ea
va ∩∩∩∩ + −
��
ib x′′d e′′b
�� xg
�
vb ∩∩∩∩ + − ∩∩∩∩
��
ic x′′d e′′c
��
�
vc ∩∩∩∩ + −
��

Figure 4: Equivalent Network Model

10 Restatement Of The Model

The synchronous machine model which uses the three phase currents as state variables may now
be stated in the form of a set of differential and algebraic equations:

dψkd
= −ω0 rkd ikd (129)
dt
dψkq
= −ω0 rkq ikq (130)
dt
dψf
= −ω0 rf if (131)
dt
dδ
= ω − ω0 (132)
dt

17
 dω ω0  
= Tm + e′′q iq + ed′′ id (133)
dt 2H
where: " # " #−1 " # " # !
ikd xkd xad ψkd xad
= − id
if xad xf ψf xad
and
1 xaq
ikq = ψkq − iq
xkq xkq
(It is assumed here that the difference between subtransient reactances is small enough to be
neglected.)
The network interface equations are, from the network to the machine:

2π 2π
id = ia cos(θ) + ib cos(θ − ) + ic cos(θ + ) (134)
3 3
2π 2π
iq = −ia sin(θ) − ib sin(θ − ) − ic sin(θ + ) (135)
3 3
and, in the reverse direction, from the machine to the network:
ω ′′
e′′a = − (e sin(θ) − ed′′ cos(θ))
ω0 q
1 de′′q 1 de′′
+ cos(θ) + sin(θ) d (136)
ω0 dt ω0 dt
ω ′′ 2π 2π
e′′b = − (eq sin(θ − ) − e′′d cos(θ − ))
ω0 3 3
1 2π de′′q 1 2π de′′d
+ cos(θ − ) + sin(θ − ) (137)
ω0 3 dt ω0 3 dt
ω 2π 2π
e′′c = − (e′′q sin(θ + ) − e′′d cos(θ + ))
ω0 3 3
1 2π de′′q 1 2π de′′d
+ cos(θ + ) + sin(θ + ) (138)
ω0 3 dt ω0 3 dt
And, of course,

θ = ω0 t + δ (139)
e′′q = ψd′′ (140)
e′′d = −ψq′′ (141)
xad xkdl ψf + xad xf l ψkd
ψd′′ = (142)
xad xkdl + xad xf l + xkdl xf l
xaq
ψq′′ = ψkq (143)
xaq + xkql

18
11 Network Constraints
This model may be embedded in a number of networks. Different configurations will result in
different constraints on currents. Consider, for example, the situation in which all of the terminal
voltages are constrained, but perhaps by unbalanced (not entirely positive sequence) sources. In
that case, the differential equations for the three phase currents would be:
x′′d dia x′′d + 2xg  xg
= (va − e′′a ) − (vb − eb′′ ) + (vc − ec′′ ) ′′

′′ (144)
ω0 dt xd + 3xg xd + 3xg
x′′d dib ′′
x + 2xg  xg
= (vb − e′′b ) d′′ − (va − ea′′ ) + (vc − ec′′ ) ′′

(145)
ω0 dt xd + 3xg xd + 3xg
x′′d dic ′′
x + 2xg xg
= (vc − e′′c ) ′′d − (vb − eb′′ ) + (va − ea′′ ) ′′
 
(146)
ω0 dt xd + 3xg xd + 3xg

12 Example: Line-Line Fault

We are not, however, constrained to situations defined in this way. This model is suitable for
embedding into network analysis routines. It is also possible to handle many different situations
directly. Consider, for example, the unbalanced fault represented by the network shown in Figure
5. This shows a line-line fault situation, with one phase still connected to the network.

′′
i�
a ra x′′d ��
ea
va ∧∧∧ ∩∩∩∩ + −
∨∨
��
ib ra x′′d e′′b
�� xg
�
∧∧∧
∨∨
∩∩∩∩ + − ∩∩∩∩
��
x′′d e′′c
��
ra
∧∧∧
∨∨
∩∩∩∩ + −
��

Figure 5: Line-Line Fault Network Model

In this situation, we have only two currents to worry about, and their differential equations
would be:
dib ω0 ′′
= (e − e′′b − 2ra ib ) (147)
dt 2x′′d c
dia ω0
= (va − e′′a − ra ia ) (148)
dt x′′d + xg
and, of course, ic = −ib .
Note that here we have included the effects of armature resistance, ignored in the previous
section but obviously important if the results are to be believed.

19
13 Permanent Magnet Machines
Permanent Magnet machines are one state variable simpler than their wound-field counterparts.
They may be accurately viewed as having constant field current. Assuming that we can define the
internal (field) flux as:
ψ0 = xad if 0 (149)

13.1 Model: Voltage Driven Machine

We have a reasonably simple expression for the rotor currents, in the case of a voltage driven
machine:
" # " #−1 " #
id xd xad ψd − ψ0
= (150)
ikd xad xkd ψkd − ψ0
" # " #−1 " #
iq xq xaq ψq
= (151)
ikq xaq xkq ψkq

The simulation model then has six states:

dψd
= ω0 vd + ωψq − ω0 ra id (152)
dt
dψq
= ω0 vq − ωψd − ω0 ra iq (153)
dt
dψkd
= −ω0 rkd ikd (154)
dt
dψkq
= −ω0 rkq ikq (155)
dt
dω ω0
= (ψd iq − ψq id + Tm ) (156)
dt 2H
dδ
= ω − ω0 (157)
dt

13.2 Curent-Driven Machine Model

In the case of a current-driven machine, rotor currents required in the simulation are:
1
ikd = (ψkd − xad id − ψ0 ) (158)
xkd
1
ikq = (ψkq − xaq iq ) (159)
xkq

Here, the “flux behind subtransient reactance” is, on the direct axis:
xkdl ψ0 + xad ψkd
ψd′′ = (160)
xad + xkdl
and the subtransient reactance is:
x′′d = xal + xad ||xkdl (161)

20
 On the quadrature axis,
xad ψkq
ψq′′ = (162)
xad + xkql
and
x′′q = xal + xaq ||xkql (163)
In this case there are only four state equations:
dψkd
= −ω0 rkd ikd (164)
dt
dψkq
= −ω0 rkq ikq (165)
dt
dω ω0  ′′ 
= eq iq + ed′′ id + Tm (166)
dt 2H
dδ
= ω − ω0 (167)
dt
The interconnections to and from the network are the same as in the case of a wound-field
machine: in the “forward” direction, from network to machine:
2π 2π
id = ia cos(θ) + ib cos(θ − ) + ic cos(θ + ) (168)
3 3
2π 2π
iq = −ia sin(θ) − ib sin(θ − ) − ic sin(θ + ) (169)
3 3
and, in the reverse direction, from the machine to the network:
ω ′′
e′′a = − (e sin(θ) − ed′′ cos(θ))
ω0 q
1 de′′q 1 de′′
+ cos(θ) + sin(θ) d (170)
ω0 dt ω0 dt
ω 2π 2π
e′′b = − (e′′q sin(θ − ) − e′′d cos(θ − ))
ω0 3 3
1 2π de′′q 1 2π de′′d
+ cos(θ − ) + sin(θ − ) (171)
ω0 3 dt ω0 3 dt
ω 2π 2π
e′′c = − (e′′q sin(θ + ) − e′′d cos(θ + ))
ω0 3 3
1 2π de′′q 1 2π de′′d
+ cos(θ + ) + sin(θ + ) (172)
ω0 3 dt ω0 3 dt

13.3 PM Machines with no damper

PM machines without much rotor conductivity may often behave as if they have no damper winding
at all. In this case the model simplifies even further. Armature currents are:
1
id = (ψd − ψ0 ) (173)
xd
1
iq = ψq (174)
xq

21
 The state equations are:
dψd
= ω0 vd + ωψq − ω0 ra id (175)
dt
dψq
= ω0 vq − ωψd − ω0 ra iq (176)
dt
dω ω0
= (ψd iq − ψq id + Tm ) (177)
dt 2H
dδ
= ω − ω0 (178)
dt

13.4 Current Driven PM Machines with no damper

In the case of no damper the machine becomes quite simple. There is no “internal flux” on the
quadrature axis. Further, there are no time derivatives of the internal flux on the d- axis. The only
machine state equations are mechanical:
dω ω0
= (ψ0 iq + Tm ) (179)
dt 2H
dδ
= ω − ω0 (180)
dt
The “forward” network interface is as before:
2π 2π
id = ia cos(θ) + ib cos(θ − ) + ic cos(θ + ) (181)
3 3
2π 2π
iq = −ia sin(θ) − ib sin(θ − ) − ic sin(θ + ) (182)
3 3
and, in the reverse direction, from the machine to the network, things are a bit simpler than before:
ω
e′′a = − ψ0 sin(θ) (183)
ω0
ω 2π
e′′b = − ψ0 sin(θ − ) (184)
ω0 3
ω 2π
e′′c = − ψ0 sin(θ + ) (185)
ω0 3
(186)

22