Power System Stability-Chapter 3

Download as pdf or txt
Download as pdf or txt
You are on page 1of 84

C H A P T E R

5
Power System Stability

5.1 Introduction
The present day electric power system network forms a complete, nonlinear
dynamical system. Several controls are provided in order to perform the power
system into a proper way. During the normal operation, all these controls try
to bring the system to an operating equilibrium ensuring the balance of real
and reactive powers in the system. Following a disturbance, the balance of
real and reactive powers gets disturbed. The dynamical power system network
undergoes a transition period and may settle down to an operating equilibrium
with the help of the above controls, which may or may not be the same as
the predisturbance equilibrium point. The capability of the system to achieve
an operating equilibrium, after disturbances, depends on its inherent strength,
nature and amount of disturbances. The system becomes unstable if it is not
capable of regaining the operating equilibrium. Thus, a general definition of
stability is given as follows:
“The stability of a dynamical system is its property or ability to remain
in a state of operating equilibrium under normal operating conditions
and to regain an acceptable state of equilibrium after being subjected to
a disturbance”.
During the early part of the 20th century, the concern of the power system
engineers was to maximize the real power transfer from the remotely located
generating stations to the load centre. The problem of maintaining synchronous
operation started when two or more generators were connected in the network to
share the system power demand. The difficulty in maintaining the synchronous
operations was experienced specifically in the case of severe disturbances
such as network faults and outage of large generating plants. This was called
the transient stability problem. Several practical measures were suggested to
improve the transient stability including the fast exciters and protection system.
Although these measures helped in improving the synchronizing capability and
the transient stability limit, a few of these resulted into deterioration in system
damping. With poor damping, the system becomes oscillatory unstable, called
217
218 Electrical Power Systems: Analysis, Security and Deregulation

small signal stability problem. Stabilizing controls such as power system


stabilizers have been used to improve the system small signal stability, facilitating
the transmission system to be utilized close to their maximum power transfer
capability. The stressed operation of the power system, due to increased real
power transfer capability, makes the transfer of reactive power difficult, giving
rise to a new stability phenomenon known as voltage stability.

5.2 Classification of Power System Stability


Although the power system stability is one single phenomenon, it has been
classified into various types for the ease of analysis and identifying the factors
affecting the stability and hence planning the control actions for its enhancement.
The classification of power system stability is shown in Figure 5.1. It is broadly
classified into two types: 1. angular stability and 2. voltage stability.

Figure 5.1 Classification of power system.

The (voltage) angle stability is the property of the interconnected power


system network to maintain synchronous operation of various generating
plants subjected to a disturbance. It is mainly concerned with maintaining real
power requirement in the system. The stability problem involves the study of
electromechanical oscillations inherent in the power systems. Disturbances in
the system can be large or small, gradual or sudden. Depending on the nature
of disturbance, the angle stability can be further classified as transient stability
or small signal stability.
Power System Stability 219
Transient stability
Transient stability refers to the large disturbances in the system. It can be
defined as the “ability of the system and its generating units to remain in
synchronism following a large (severe) and sudden disturbance”. Faults
in the transmission system, sudden change of bulk load, loss of operating
units, line switching are the examples of large disturbances. In general, the
postdisturbance operating equilibrium is different from the predisturbance
equilibrium point in case of such disturbances. Since the severe disturbances
involve a large deviation in rotor angles, nonlinear dynamical model of the
system is considered for the transient stability studies.

Small signal stability


Small signal stability refers to the “ability of the system to maintain
synchronism under small and sudden disturbances”. Such disturbances occur
continuously in the system due to small variation in loads and generations.
Small signal instability can be either due to insufficient synchronizing torque
resulting into monotonous increase in the rotor angles or due to insufficient
damping torque resulting into undamped angle oscillations in the system.
The stability depends on several factors such as initial operating point,
transmission system strength, generator excitation and other controls in the
system. Since the disturbances are small, linearized dynamical model of the
system can be used for analysis. The small signal stability can be further
divided into the following types.
(i) Local modes or machine-system mode, due to swinging of a generating
plant unit with the rest of the system. The frequency of oscillation
may range from 0.7 to 2 Hz.
(ii) Interarea mode, due to swinging of a group of coherent generating
units with other group(s) of coherent units. In general, these groups
are interconnected with weak tie lines. The frequency of oscillation
may range from 0.1 to 0.7 Hz.
(iii) Control mode, due to the poorly tuned controllers in the system such
as exciters, speed governors, HVDC converters, SVC, etc.
(iv) Torsional mode, due to interaction of mechanical oscillation of turbine
generator shaft system with the oscillations in the electrical circuit,
involving various controls and series compensated lines.

Voltage stability
The voltage stability also called load stability refers to the “ability of the
system to maintain load bus voltages within acceptable limit, following
some disturbance or change in power demand”.
The voltage stability can be further classified as follows.
(i) Large disturbance voltage stability, which is the ability of the system
to regain voltages at all the buses within the acceptable steady state
220 Electrical Power Systems: Analysis, Security and Deregulation

levels, when subjected to a large disturbance. The study ranges from


a few seconds to tens of minutes and requires nonlinear dynamic
simulation of load characteristics, and various controls including load
tap changer dynamics and generator current limiters.
(ii) Small disturbance voltage stability considers small disturbances such
as incremental change in loads. Static analysis is utilized to identify
voltage instability conditions, various contributing factors and the
stability margin.
An alternate classification of power system transient stability study is
based on the time frame considered and the time response of the concerned
phenomenon. It is broadly classified into three types.
(i) Short-term for a time period between 0 and 10 seconds
(ii) Mid-term for 10 seconds to a few minutes
(iii) Long-term for a few minutes to tens of minutes.
The modelling details of various power system components depend on
the time frame of the phenomenon being studied. For example, the network
transient is ignored in mid-term and long-term stability studies. Slow dynamics
such as boiler, and on-load tap changer dynamics is required to be considered
in the long-term stability study.

5.3 Transient Stability


5.3.1 Introduction
The power system transient stability has been defined as the ability of the
synchronous generators in an interconnected network to remain in synchronism
after being subjected to a large or severe disturbance. During a steady state
operation, the generators run at a constant synchronous speed, with rotor
acceleration being zero, maintaining a balance between the mechanical input
power (Pm) from the turbine and the electrical power (Pe ) output to the
electrical network. A disturbance in the electrical network causes the (real
power) output of generators to change creating an imbalance with the
mechanical power input. Since the mechanical power input from turbine cannot
change instantaneously, the change in electrical power requirement is met
initially from the stored kinetic energy of the rotors causing the rotors to
accelerate or decelerate and also changing the rotor angle position. Severe
disturbances may cause a large excursion in the rotor angles. The transient
stability requires
(i) the calculation of the electrical power output of generators during
prefault, fault and postfault conditions. This involves solution of
network power flow equations, which are nonlinear algebraic in
nature.
Power System Stability 221
(ii) solving the rotor dynamic equations to study the variation of rotor
angle with time. The generator mechanical dynamics is described by
the nonlinear differential equation known as swing equation. This
requires some numerical integration method for solution.
The system is said to be stable during the postdisturbance period, if the
rotors of all the machines achieve constant synchronous speed. To simplify
the transient stability analysis, the classical approach has been utilized which
is also discussed in this chapter.

5.3.2 Assumptions of Transient Stability Analysis


The classical technique involves the following main assumptions.
1. Mechanical input to the generator remains constant (i.e. speed governing
system action is neglected).
2. Machine damping and AVR action are neglected. Synchronous machine
is modelled as constant voltage source behind the transient reactance.
3. Network transients are neglected. Thus static model of network can be
used.
4. Loads are represented as constant impedance/admittance type.
5. Mechanical angle of each machine rotor coincides with the electrical
phase of voltage behind transient reactance.
In addition, the transmission line resistances and saliency of the synchronous
machine can also be neglected, which gives conservative results. The static
network equations and the machine dynamics can be formulated as given below.

5.4 Power Angle Equation of a Two


Machine System
Consider a very simple system as shown in Figure 5.2.

Figure 5.2 Two machine system.

It consists of a synchronous generator supplying power to a synchronous


motor through a transmission line.
To draw the reactance diagram, we know that any synchronous machine
is represented by a constant voltage source in series with a reactance X.
Depending upon the condition under study, the reactance may be subtransient
reactance Xd, transient reactance Xd or steady state synchronous reactance Xd.
Thus the above one line diagram can be drawn in terms of reactance diagram
as shown in Figure 5.3.
222 Electrical Power Systems: Analysis, Security and Deregulation

Figure 5.3 Reactance diagram.

In Figure 5.3, the generator is represented by EG in series with XG,


the motor is represented by EM in series with XM and
the transmission line with leakage reactance is given by XL.
The total reactance between the machines is given by
X = XG + XL + XM (5.1)
The internal voltages EG and EM are generated by the flux produced by the
field windings of the machines, hence their phase difference is the same as
the electrical angle between the machine rotors.
The vector diagram is shown in Figure 5.4.

Figure 5.4 Vector diagram.

From the vector diagram,


EG = EM + jX Ir (5.2)
The equation for current is

EG - EM
I = (5.3)
jX
Since the resistances of the machines and the transmission lines are neglected,
the power output of the generator is also the power input to the motor and
is given by
*
P = Real part of ( EG ¥ I )

= Re ; EG c G mE
* E - EM
(5.4)
jX
Power System Stability 223

Let EM = EM0, EG = EGd,


*
 EG = EG–d,
Upon substitution,
È Ê E –d - EM –0 ˆ ˘
P = Re Í EG –-d Á G ˜¯ ˙
Î Ë X –90∞ ˚
È E2 E E ˘
= Re Í G –-90∞ - G M –-90∞ - d ˙ (5.5)
Î X X ˚
 Real power P is given by
EG EM
P= -
cos ( -90∞ - d )
X
EG EM
= - cos (90∞ + d )
X
E E
P = G M sin d (5.6)
X
The above equation shows that the power transmitted from the generator to
the motor varies with the sine of the displacement angle d between the two
rotors. Hence this equation is called the power angle equation and if the curve d
against P is plotted, it is called the power angle curve as shown in Figure 5.5.

Figure 5.5 Power angle curve.

The maximum power occurs when sin d = 1, i.e. d = 90°.


EG EM
 Pmax = is the steady state stability limit.
X
When the slope dP/dd is positive (–90°  d  90°), it means that an increase
in displacement angle results in an increase in transmitted power and hence
the system will be stable. If dP/dd is negative, it indicates that the system is
unstable.
224 Electrical Power Systems: Analysis, Security and Deregulation

EXAMPLE 5.1 Two synchronous machines of equal rating have internal


voltages of 1.1 + j0.5 and 0.8 – j0.4 per unit voltages respectively. The machines
are connected by a line of 50 km length having only reactance and the second
machine receives power of 0.9 per unit. Determine the reactance of the line
per km length. Assume that there is no internal reactance for simplification.
Solution:
Given
EG = 1.1 + j0.5 = 1.2124.4°
EM = 0.8 – j0.4 = 0.89–26.6°
P = 0.9, length of transmission line = 50 km
XG = XM = 0
d = 24.4° – (–26.6°) = 51°
We know that the power angle equation is
EG EM
P= sin d
X
1.21 ¥ 0.89
0.9 = sin 51∞
X
0.8087
X= = 0.8986 p.u.
0.9
Since XG = XM = 0, X denotes the reactance of the transmission line.
X = 0.8986 p.u
0.8986
X per km = = 0.0186 p.u.
50

5.5 Power Angle Equation of a Salient


Pole Machine
Consider a salient pole machine represented by means of a one line diagram
as shown in Figure 5.6.

Figure 5.6 Representation of a salient pole machine.

Let Ef be the excitation emf per phase


Vt be the terminal voltage per phase
Vb be the voltage of the infinite bus
Power System Stability 225
Xd be the direct axis synchronous reactance
Xq be the quadrature axis synchronous reactance
Xext be the reactance between generator and infinite bus
The electrical power output of a salient pole generator is given by,
E f Vb Vb2 ( X d - X q )
Pe = sin d + sin 2d (5.7)
Xd 2 Xd Xq
( Excitation power ) (Reluctance power)
The excitation power is the same as the power angle equation of a simple
two machine system. The reluctance power varies as sin 2d with the maximum
value of d = 45°. This term is independent of field excitation and would be
present even if the field is unexcited. The reluctance component is of the order
of 10 to 20 per cent of the excitation component. The reluctance component
of power is usually neglected in the steady state stability studies.
For a non-salient pole machine, Xd = Xq, we get the original power angle
equation of a simple two machine system.

5.6 Swing Equation


The behaviour of a synchronous machine during the transient period is described
by the swing equation.
We know that the torque exerted on a rotating body is given by the product
of moment of inertia J (kg · m2) and angular acceleration a (rad/s2), that is,

d 2q
Ta = Ja = J (5.8)
dt 2
where q is the angular position of the rotor in radians at any instant of time,
and t is the time in seconds.
It is convenient to measure q with respect to a reference axis that is rotating
at the synchronous speed. If d is the angular displacement of the rotor in
electrical degrees from the synchronously rotating reference axis and ws is the
synchronous speed in electrical radians, then q can be expressed as the sum
of: (1) time varying angle wst on the rotating reference axis, and (2) the torque
angle d of the rotor with respect to the rotating reference axis. In other words,
q = wst + d   electrical radians (5.9)
Differentiating with respect to t, we get
dq dd
= ws + (5.10)
dt dt
Differentiating once again with respect to t, we get

d 2q d 2d
= (5.11)
dt 2 dt 2
226 Electrical Power Systems: Analysis, Security and Deregulation

d 2q
 Angular acceleration of rotor, a =
dt 2
d 2d
a= electrical radians (5.12)
dt 2
In a synchronous generator, the accelerating torque Ta is equal to the difference
of input shaft torque Tm and the output electromagnetic torque Te.

 Ta = Tm – Te

d 2q d 2d
J◊ =J◊ (5.13)
dt 2 dt 2
d 2d
 Ta = Tm – Te = J ◊
dt 2
Multiplying both sides by angular velocity, w, we get

d 2d
wTa = wTm - wTe = J w (5.14)
dt 2
If Pa, Pm and Pe denote the accelerating power, mechanical power input and
electrical power output respectively, we get

d 2d
Pa = Pm - Pe = J w (5.15)
dt 2
Also, since angular momentum M = Jw, therefore Eq. (5.15) can be written as

d 2d
Pa = Pm - Pe = M (5.16)
dt 2
Equation (5.16) is called the swing equation. It is a nonlinear differential
equation of the second order.

5.6.1 Swing Curves


The solution of swing equation gives the variation of d (in electrical radians)
with respect to time (in seconds). The graph when plotted is called a swing
curve as shown in Figure 5.7. It provides information regarding stability.
They show the tendency of d to oscillate and/or increase beyond the point of
return. If d increases continuously with time, the system is unstable. While
if d starts decreasing after reaching a maximum value, it is inferred that the
system will remain stable.
Power System Stability 227

Figure 5.7 Swing curves.

5.6.2 Constants Used in Stability Analysis


Inertia constant, M
d 2d
From the swing equation M = Pa
dt 2
If power is in W, d is in rad and t is in s, then M is in W/rad/s2 or
W· s2/rad. Since 1 J = 1 W· s, M can have its unit as J · s/rad. If power is in
MW, then M is MJ · s/rad.
If d is specified in electrical degrees, then M has unit as MJ · s/electrical
degree.
If power is in pu, then M has pu power s2/electrical degree. The value of
M will be referred to as per unit value of M.
In general, constant M may be defined as the power in MW required to
produce unit angular acceleration.
Kinetic energy, N
The kinetic energy of the rotor at synchronous speed denoted by N is given by
1
N = Mw s
2
where ws = 2p ns ; mechanical rad/s
= 2p f ; electrical rad/s
= 360 f ; electrical degree/s
where ns is the speed in rps and f is the frequency in Hz.
Normally generators of the same MVA ratings may have different values
of kinetic energy and momentum. To express them in a common way, we use
a constant H (also called inertia constant).
Inertia constant, H
It is defined as the ratio of stored kinetic energy to volt ampere rating of
machine.
kinetic energy
H= MJ/MVA
MVA rating
228 Electrical Power Systems: Analysis, Security and Deregulation

N
= where S is the MVA rating
S
 N = SH

Equivalent H constant
Consider a system in which ‘n’ number of generators are connected in parallel
to the same bus bar.
Let S1, S2, S3, …, Sn be the MVA rating of individual machines
H1, H2, H3, …, Hn be the inertia constants of individual machines
N1, N2, N3, …, Nn be the kinetic energy stored in individual machines
Se be the MVA rating of equivalent machine
He be the inertia constant of equivalent machine
Ne be the kinetic energy stored in equivalent machine and Sb be the
base MVA.
The energy stored by the equivalent machine is given by the sum of energies
stored by individual machines.
N = N1 + N2 +  + Nn
Se He = S1H1 + S2H2 +  + SnHn
where Se = S1 + S2 +  + Sn
If the base MVA, Sb, is equal to the combined MVA rating of individual
machines Se , i.e. Sb = Se, we get
ÊS ˆ ÊS ˆ ÊS ˆ
H e = H1 Á 1 ˜ + H 2 Á 2 ˜ +  + H n Á n ˜
Ë b¯
S Ë b¯
S Ë Sb ¯

If the machines are identical, we have
S1 = S2 =  = Sn = S
H 1 = H 2 =  = Hn = H
then
Ê HS ˆ
He = n ¥ ÁË ˜
nS ¯
With identical machines
Sb = Se = n  S
Substituting Sb, we get
Ê HS ˆ
He = n ¥ Á
Ë nS ˜¯
He = H
Thus the equivalent H constant of several identical machines operating in
parallel is the same as that of any one of the machines.
Power System Stability 229
Equivalent M constant of two machines
Two synchronous machines connected by a reactance can be replaced by one
equivalent machine connected through a reactance to an infinite bus as follows.
The swing equation of machine 1 is given by
d 2d1
M1 = Pm1 - Pe1 (5.17)
dt 2
The swing equation of machine 2 is given by
d 2d 2
M2 = Pm 2 - Pe 2 (5.18)
dt 2
From Eq. (5.17)
d 2d1 Pm1 - Pe1
2
= (5.19)
dt M1
From Eq. (5.18)
d 2d 2 Pm 2 - Pe 2
= (5.20)
dt 2
M2
Subtracting Eq. (5.20) from Eq. (5.19), we get

d 2d1 d 2d 2 Pm1 - Pe1 Pm 2 - Pe 2


- = -
2 2
dt dt M1 M2
We can write
d 2 (d1 - d 2 ) M 2 ( Pm1 - Pe1 ) - M1 ( Pm 2 - Pe 2 )
= (5.21)
2
dt M 1M 2
If d is the relative angle between the rotors of the two machines, then d =
d1 – d2. We can write Eq. (5.21) as
d 2d M 2 Pm1 - M1 Pm 2 M 2 Pe1 - M1 Pe 2
= - (5.22)
2
dt M 1M 2 M 1M 2
Multiplying both sides of the above equation by M1M2 /(M1 + M2), we get

M1M 2 d 2d M1M 2 Ê M 2 Pm1 - M1 Pm 2 M 2 Pe1 - M1 Pe 2 ˆ


2 =
-
M1 + M 2 dt M1 + M 2 ÁË M 1M 2 M 1M 2 ˜¯

M1M 2 d 2d M P - M1 Pm 2 M 2 Pe1 - M1 Pe 2
= 2 m1 - (5.23)
M1 + M 2 dt 2 M1 + M 2 M1 + M 2
The swing equation of an equivalent machine is given by
d 2d
= Pm¢ - Pe¢
M¢ (5.24)
dt 2
From Eq. (5.24) we can say that the equivalent values of M, Pm and Pe are
given by
230 Electrical Power Systems: Analysis, Security and Deregulation

M 1M 2
M =
M1 + M 2
M 2 Pm1 - M1 Pm 2
Pm =
M1 + M 2
M 2 Pe1 - M1 Pe 2
and Pe =
M1 + M 2
Relationship between inertia constant M and inertia constant H
We have
2N 2N N
M= = =
w s 360 f 180 f
SH
= MJ · s/electrical degree   (since N = SH)
180 f
If the angle is expressed in radians
SH
M= MJ · s/electrical radian
pf
EXAMPLE 5.2 The moment of inertia of a 4 pole, 100 MVA, 11 kV,
3-f, 0.8 power factor, 50 Hz turbo alternator is 10000 kg · m2. Calculate H
and M.
Solution:
J = 10000 kg · m2
120 f 120 ¥ 50
Ns = = = 1500 rpm
p 4
N s 1500
ns = = = 25 rps
60 60
ws = 2p ns = 50p
1 1
N= J w s2 = ¥ 10000 ¥ (50p ) 2 = 123.37 MJ
2 2
N 123.37
H= = = 1.2337 MJ/MVA
S 100
SH 100 ¥ 1.2337
M= = = 0.0137 MJ ◊ s /electrical degree
180 f 180 ¥ 50
EXAMPLE 5.3 A 50 Hz, 4 pole, turbo alternator rated 100 MVA, 11 kV
has an inertia constant of 8 MJ/MVA. Determine
1. the energy stored in the rotor at synchronous speed.
2. find the rotor acceleration if the mechanical input is suddenly raised
to 80 MW for an electric load of 50 MW. (Neglect mechanical and
electrical losses).
Power System Stability 231
Solution:
(i) H = 8 MJ/MVA; S = 100 MVA
We know that
N = HS = 800 MJ
d 2q
(ii) Swing equation is M = Pa = Pm - Pe
dt 2
Here for alternator, Pm = 80 MW
Pe = 50 MW
 Pa = 30 MW

SH N 800
Also M = = = = 0.0889 MJ◊ s
180 f 180 f 180 ¥ 50
d 2d Pa 30
 Acceleration= = = = 337.5 electrical degree/s
dt 2
M 0.0889
EXAMPLE 5.4 A 50 Hz, 4 pole turbo generator rated 20 MVA, 11 kV has
an inertia constant of H = 9 kW· s/kVA. Find the kinetic energy stored in the
rotor at synchronous speed. Find the acceleration, if the input less the rotational
losses is 26800 HP and the electrical power developed is 16 MW.
Solution:
S = 20 MVA, 11 kV,
H = 9 kW· s/kVA = 9 kJ/kVA
H = 9 MJ/MVA
Kinetic energy = N = HS = 9  20 = 180 MJ
Pa = Pm – Pe
Pm = 26800  746 = 19992800 W = 19.99 MW
Pa = 19.99 – 16 = 3.99 MW

N 180
M= = = 0.02 MW ◊ s 2 / electrical degree
180 f 180 ¥ 50
We know that
d 2d d 2d Pa 3.99
M 2
= Pa fi 2
= =
dt dt M 0.02
d d 2
Acceleration 2
= 199.5 electrical degree/s 2
dt
EXAMPLE 5.5 A power station with 4 generators each 80 MVA, 8 MJ/MVA
is in proximity with another power station having 3 generators each 200 MVA,
3.5 MJ/MVA. Determine the inertia constant of a single equivalent machine
for use in stability studies. Assume a base value of 100 MVA.
232 Electrical Power Systems: Analysis, Security and Deregulation

Solution:
4 generators are each of 80 MVA, 8 MJ/MVA
3 generators are each of 200 MVA, 3.5 MJ/MVA
We know that
HeSe = H1S1 + H2S2 + 
4 3
= Âi =1
H1S1 + ÂH S
i =1
2 2

= (4  80  8) + (3  200  3.5)
= 2560 + 2100 = 4660 MJ
4660 4660
He = = = 46.6 MJ/MVA
Se 100
EXAMPLE 5.6 Two turbo alternators specified below are interconnected
using a short line:
Machine 1 : 4 poles, 50 Hz, 125 MVA, 0.8 lag, 25000 kg · m2
Machine 2 : 4 poles, 50 Hz, 150 MVA, 0.9 lag, 20000 kg · m2
Determine the inertia constant of the single equivalent machine on a base of
150 MVA.
Solution:
Machine 1
N1
M1 =
180 f
1
N1 = J w s2
2
120 ¥ 50
Ns = = 1500 rpm
4
1500
ns = = 25 rps
60
1
 N1 = ¥ 25000 ¥ w s2
2
ws = 2p ns = 157.0796 rad/s = 308.425 MJ
N1
M1 = = 0.03426 MJ ◊ s/electrical degree
180 f
M1 in p.u. on a base of 150 MVA = 2.2846  10–4 p.u.
Machine 2
ns = 25 rps
1 2 1 2
N2 = J w s = ¥ 20000 ¥ (157.0796) = 246.74
2 2
Power System Stability 233

N2
M2 = = 0.274155 MJ ◊ s/electrical degree
180 f
M2 in p.u on a base of 150 MVA = 1.8277  10–4 p.u.
M 1M 2
M = = 1.01538 ¥ 10-4 p.u.
M1 + M 2

EXAMPLE 5.7 A generator A is rated at 50 Hz, 60 MW, 75 MVA,


1500 rpm and has an inertia constant H = 7 MJ/MVA. The corresponding
data for another generator B is 50 Hz, 120 MW, 133.3 MVA, 3000 rpm, and
4 MJ/MVA.
(a) If these two generators operate in parallel in a power station, calculate
H for the equivalent generator on a base of 100 MVA.
(b) If the power station is connected to another power station which has
two of each type of generator, calculate H for the equivalent generator
connected to an infinite bus bar.
Solution:
(a) We know that
Se He = H1S1 + H2S2
He Se = (7  75) + (4  133.3)
(7 ¥ 75) + (4 ¥ 133.3)
He = = 10.582 MJ/MVA
100
(b) Another power station has two of each type of generator. So the
equivalent He becomes twice the original He1.
He2 = 2  10.582 = 21.164 MJ/MVA
When they are connected in parallel, the equivalent machine
H e1 ◊ H e 2
He = = 7.055 MJ/MVA
H e1 + H e 2

EXAMPLE 5.8 Two turbo alternators given below are interconnected using
a short line.
Machine 1 : 4 poles, 50 Hz, 75 MVA, 0.8 lag, 30000 kg · m2
Machine 2 : 2 poles, 50 Hz, 100 MVA, 0.85 lag, and 10000 kg · m2
Determine the inertia constant of the single equivalent machine on a base of
200 MVA.
Solution:
Machine 1
120 f
ws = 2p ns = 2p ¥ = 157.079
4 ¥ 60
234 Electrical Power Systems: Analysis, Security and Deregulation

1 1
Kinetic energy = J w s2 = ¥ 30000 ¥ (157.079)2
2 2
N1 = 370
N1
H1 = = 4.935 MJ/MVA
S1
Machine 2
ws = 314.159
1 1
N2 = J w s2 = ¥ 10000 ¥ (314.159) 2 = 493.48
2 2
N 2 493.48
H2 = = = 4.9348 MJ/MVA
S2 100
Equivalent inertia constant
75
¥ 4.935
SH1
M1 = = 200 = 2.056 ¥ 10-4 p.u.
180 f 180 ¥ f
SH 2
M2 = = 2.738 ¥ 10-4 p.u.
180 f
M 1M 2
M= = 1.174 ¥ 10-4 p.u.
M1 + M 2

5.6.3 Determination of Change in Rotor Angle


When the Machine is Loaded
The swing equation is given by
d 2d
M = Pa
dt 2
d 2d Pa
 =
dt 2
M
where Pa /M is a constant.
dd
Multiplying both sides by 2
dt
d d d 2d Pa dd
= ◊2 2
dt dt 2 M dt
2Pa d d
= ◊
M dt
2
d È dd ˘ P dd
Í ˙ = 2 a
dt Î dt ˚ M dt
Power System Stability 235
Upon integration,
2
Ê dd ˆ P
Ú dÁ ˜ = 2 a
Ë dt ¯ M Ú dd
2
Ê dd ˆ Pa
ÁË ˜¯ = 2 ◊ d
dt M
dd 2 Pa 1/ 2
= d
dt M
dd 2Pa
1/ 2 = ◊ dt (5.25)
d M
Integrating once again,

Úd Ú
-1/ 2 2Pa
dd = ◊ dt
M
d -1/ 2 +1 2Pa
= ◊t
1 M
- +1
2
d 1/ 2 2Pa
= ◊t (5.26)
1 M
2
1 2Pa Pa
d 1/ 2 = 2 M ◊ t = 2M
◊t

Substituting Eq. (5.26) in Eq. (5.25), we get

dd 2 Pa Pa 2 Pa2
= ◊ ◊t = ◊t
dt M 2M 2M 2
dd ÊP ˆ
= Á a ˜ t is the change in rotor angle at anytime
dt ËM¯

EXAMPLE 5.9 The rotor of an alternator is subjected to an acceleration of


15 electrical rad/s2. If this acceleration exists constantly, compute the change
in rotor angle and the speed in rpm at the end of 5 cycles. H = 5 MJ/MVA.
Frequency = 50 Hz. Initially the machine is running at a normal speed without
any acceleration. Number of poles = 2.
Solution:
d 2d
Given the acceleration 2 = 15 electrical rad/s2
dt
d 2d Pa
From swing equation, = = 15 electrical rad/s2
dt 2
M
236 Electrical Power Systems: Analysis, Security and Deregulation

d d Pa
Change in rotor angle = = ◊t
dt M
Time period for five cycles = no. of cycles  time period for the given frequency
1
= 5 ¥ = 0.1 s
50
dd
 Change in rotor angle = = 15  0.1 = 1.5 electrical rad/s
dt
dd
= 1.5 electrical rad/s
dt
Since pole pair p = 1, we have qm = qe
dd
 = 1.5 mechanical rad/s = 2p ns
dt
1.5
ns = = 0.2387 rps
2p
 Ns = 0.2387  60 = 14.32 rpm
Since the machine is subjected to constant acceleration, the speed will increase
by 14.32 rpm.
 New speed = synchronous speed + Ns
120 f
Synchronous speed = = 3000 rpm
p
 New speed = 3014.32 rpm
EXAMPLE 5.10 A 200 MVA, 11 kV, 50 Hz, 4 pole turbo alternator has an
inertia constant of 6 MJ/MVA.
(i) Determine the stored energy in the rotor at synchronous speed.
(ii) The machine is operating at a load of 120 MW when the load suddenly
increases to 160 MW. Determine the rotor retardation. Neglect losses.
(iii) The retardation calculated above is maintained for 5 cycles. Determine
the change in power angle and the rotor speed in rpm at the end of
this period.
Solution:
H = 6 MJ/MVA, S = 200 MVA
(i) N = HS = 6  200 = 1200 MJ
(ii) Pa = 40 (retardation)
N 1200
M= = = 0.1333 MJ ◊ s/ electrical degree.
180 f 180 ¥ 50
d 2d Pa 40
= = = 300 electrical degree / s 2
dt 2
M 0.1333
Power System Stability 237
We know that change in rotor angle
d d Pa
= ◊t
dt M
where t is the time period for acceleration or retardation.
1
(iii) Here for five cycles t = 5 ¥ = 0.1 s
50
dd
 = 300  0.1 = 30 electrical degree/s
dt
p
= 30 ¥ electrical rad/s = 0.5236 electrical rad/s
180
0.5236
= mechanical rad/s
p
Pole pairs, p = 2
dd 0.5236
 = = 0.2618 mechancial rad/s = 2pns
dt 2
0.2618
ns = = 0.04167 rps
2p
Ns = ns  60 = 2.5 rpm
Since the retardation is for 5 cycles,
New speed = synchronous speed – Ns
 New speed = 1500 – 2.5 = 1497.5 rpm

5.7 Equal Area Criterion


In a system to determine whether a power system is stable after a disturbance,
it is necessary to plot the swing curve. If this curve shows that the angle
between any two machines tends to increase without limit, the system is unstable.
If after all disturbances, the angle between the two machines reaches the maximum
value and thereafter decreasing, thus oscillating with constant amplitude,
it is probable although not certain that the system is stable. There is a simple
graphical method of determining whether the machines come to rest with
respect to each other. This method is known as the equal area criterion for
stability.
The principle by which stability under transient conditions is determined
without solving the swing equation is called the equal area criterion of stability.
When this method is used, it completely eliminates the need of computing
swing curves and hence saves a considerable amount of work. The derivation of
the equal area criterion is made for one machine and an infinite bus although
the method can be adapted to a two-machine system.
238 Electrical Power Systems: Analysis, Security and Deregulation

Although not applicable to a multi-machine system, the method helps


in understanding how certain factors influence the transient stability of any
system.

5.7.1 Assumptions Made in Equal Area Criterion


1. Constant input over the time interval being considered.
2. Damping effect is neglected.
3. Constant voltage behind transient reactance.
In a two-machine system, it is quite likely to be stable, if it survives the
first swing.

Validity of the assumptions


1. The assumption of constancy of the input power is justified by the
fact that the effect of the action of governors on the input will not be
appreciable up to the first swing, and thereafter, it will only tend to
restore stability.
2. The effect of damping would be to slightly reduce the amplitude of
the first swing and would therefore tend to diminish the amplitude of
any subsequent swing.
When a fault occurs in a system, the fault current tends to cause
demagnetization of the field due to armature reaction and the field current
increases to the extent required to offset this effect so as to keep up the
constant flux linkage of the field circuit. If the machine is not provided with
a voltage regulator, the field current decreases to its original value, as also
the flux linkage. As the time constant may be a few seconds, during the first
swing there will not be any appreciable decrease in the flux linkage. Even if
a voltage regulator is provided, it will take some time before the regulator and
the exciter response becomes effective. During the subsequent swings, however,
the voltage regulator will contribute substantially toward the maintenance of
stability. It is, therefore, reasonable to assume constant voltage behind transient
reactance without any significant error.

5.7.2 Equal Area Criterion Method Applied to a Machine


Swinging with Respect to an Infinite Bus
The swing equation of a synchronous machine swinging with respect to an
infinite bus is given by
d 2d
M = Pa = Pm - Pe
dt 2
d 2d Pa
=
dt 2
M
Power System Stability 239

dd
Multiplying both sides by 2
dt

d d d 2d P dd
2 =2 a
dt dt 2
M dt
d ÈÊ dd ˆ ˘
2
P dd
ÍÁ ˜ ˙ = 2 a
dt Î Ë dt ¯ ˚ M dt
ÈÊ dd ˆ 2 ˘ P
d ÍÁ ˜ ˙ = 2 a dd
Ë
Î dt ˚¯ M
Upon integration,
2 dt
Ê dd ˆ Pa
ÁË ˜¯ =
dt
d0
Ú 2 M dd
dt
2
=
M Ú P dd
d0
a

dt
dd 2
dt
=
M Ú P dd
d0
a


If dd/dt = w, the velocity of the displacement angle with respect to the infinite
bus, then
dt
dd 2

w=
dt
=
M Ú P dd
d0
a

If the machine is continuously swinging, then the above equation will be non-
zero. The stability is indicated by the zero value, i.e. if the integral vanishes.

i.e. Ú P dd
a =0

Ú (P m - Pe ) d d = 0

Ú P dd - Ú P dd = 0
m e

 Ú P dd = Ú P dd
m e

Under the steady operating conditions, Pm = Pe as given by the constant input


line. The power angle curve is shown in Figure 5.8.
Let us assume that the initial angle be d0. The corresponding value of Pe
is given by bf.
240 Electrical Power Systems: Analysis, Security and Deregulation

Figure 5.8 Equal area criterion.

For the stability criterion Ú P dd = Ú P dd


m e
which is indicated by the equality of the area of the rectangle afgd and the
area bcegf.
Area afgda = Area bcedgf b
Area abc + Area bcdgf = Area ced + Area bcdgf
Area abc = Area ced
A1 = A2 which is shown by the dashed line, i.e. the area A1 below Pm line is
equal to the area A2 above Pm line. Hence this method is called equal area
criterion.
From Figure 5.8, ‘ab’ represents Pm – Pe  accelerating power corresponding
to the initial angle d0, ‘de’ represents Pe – Pm  decelerating power and when
the accelerating power is equal to the decelerating power, the machine tends
to return to the original condition if there is a balance between the two areas
at d = dm .
The system will remain stable only when A2  A1, i.e. if accelerating area
A1 is greater than the decelerating area A2, the system will definitely become
unstable.

5.7.3 Sudden Change in Mechanical Input


Some of the conditions caused by the sudden
increase in the mechanical load on a synchronous
motor connected to an infinite bus can be predicted
by analyzing Figure 5.9. An infinite bus is
considered as a generator with infinite H constant
and with constant frequency.
The sinusoidal curve Pe is a plot of the electric
power input to the motor with all resistance Figure 5.9 Motor connected
neglected. The curve Pe is plotted from the equations to infinite bus.
Power System Stability 241

P = (|EG | | E M| / |X |) sin d and Pmax = | EG | | EM | | X | , where |EM | is the voltage of


the infinite bus, |EG | is the voltage behind transient reactance of the motor, and
X is determined from the transient reactance of the motor plus the reactance
of the transformer and line, if any, between the motor and the infinite bus.
Originally the motor is operating at synchronous speed with a torque angle
of d0, and the mechanical power output P0 is equal to the electric power input
Pe corresponding to d0. When the mechanical load is suddenly increased so
that the power output is Ps , which is greater than the electric power input
at d0, the difference in power must come from the kinetic energy stored in
the rotating system. This can be accompanied only by a decrease in speed,
which results in an increase in the torque angle d. As d increases, the electric
power received from the bus increases until Pe = Ps at point “b” in the curve.
At this point there is equilibrium of input and output torque so that acceleration
is zero, but the motor is running at less than synchronous speed so that d is
increasing. The angle d continues to increase, but after passing through point
“b” the electric power input Pe is greater than Ps, and the difference must be
stored in the system through an increase in kinetic energy accompanying an
increase in speed. Thus, between points “b” and “c” as d increases, the speed
also increases, until synchronous speed is again reached at point “c”, where the
torque angle is dm. At point “c”, Pe is still greater than Ps and speed continues
to increase, but d starts to decrease as soon as the speed of the motor exceeds
synchronous speed. The maximum value of d is dm at point “c”. As d decreases,
point “b” is reached again with the speed more than the synchronous speed,
so that d continues to decrease until point “b” is reached. The motor is again
operating at synchronous speed, and the cycle is repeated.
When the load is suddenly increased from P0 to Ps , the motor oscillates
around the equilibrium torque angle ds between d0 and dm as shown in
Figure 5.10. If damping is present, the oscillations decrease, and stable operation
results at ds. Table 5.1 shows the changes in speed, angle, electric power input,
mechanical power output, stored energy and acceleration or deceleration as
the machine oscillates. A thorough study of Table 5.1 will lead to a better
understanding of transient disturbances.

Figure 5.10 Electric power input to a motor as a function of torque angle d.


242 Electrical Power Systems: Analysis, Security and Deregulation

Table 5.1 Changing conditions in a synchronous motor swinging with respect to


an infinite bus because of a sudden increase in load
Position in Motor Torque Electric Stored Rotating
cycle speed w angle d power Pe energy system
1/2 Jw2 = undergoing
W
At point a w = ws, d = d0, Pe < Ps, w = w s, Deceleration
decreasing minimum minimum decreasing
From a w = ws, Increasing Pe < Ps, w < w s, Deceleration
towards b decreasing increasing decreasing
At point b w = ws, d = ds¸ Pe = Ps, w < w s, Acceleration
minimum increasing increasing minimum
From b w = ws, Increasing Pe > Ps, w < w s,
towards c increasing increasing decreasing
At point c w = ws, d = dm, Pe > Ps, w = w s, Acceleration
increasing maximum maximum increasing
From c w = ws, Decreasing Pe > Ps, w > w s,
towards b maximum decreasing increasing
At point b w = ws, d = ds, Pe = Ps, w > w s, Deceleration
maximum decreasing decreasing maximum
From b w = ws, Decreasing Pe < Ps, w > w s,
towards a decreasing minimum decreasing

The maximum swing of the motor to a torque angle dm can be found by


equal area criterion by equating the shaded areas A1 and A2.
The shaded area A1 is given by
ds

A1 = Ú ( P - P ) dd
d0
s e

Similarly, the shaded area A2 is given by


dm

A2 = Ú ( P - P ) dd
ds
e s

ds dm

and A1 – A2 = Ú ( P - P ) dd - Ú ( P - P ) dd
d0
s e
ds
e s

dm

A1 – A2 = Ú ( P - P ) dd
d0
s e

d
2( Ps - Pe )
Equation Ú
d0
M
d d = 0 is satisfied and dd /dt = 0 when A1 = A2. The

maximum torque angle dm is located graphically so as to make A2 equal to A1.


Power System Stability 243
5.7.4 Estimating Transient Stability Limit
Figure 5.11 shows a suddenly applied load, which is larger than that shown
in Figure 5.12. The area A2 above Ps under the curve Pe is less than A1 and
dd/dt is not zero at d = dm. Therefore d continues to increase after d = dm.
Pe again becomes less than Ps. The torque angle d continues to increase
beyond dm, and restoring forces are not encountered. The system is stable
only if an area A2 located above Ps is equal to A1. The test of equal areas is
called the equal area criterion. The maximum allowable increase in the power
suddenly taken from the motor originally supplying the power P0 is shown in
Figure 5.12. A suddenly applied load greater than that shown in Figure 5.12
would not permit the torque angle of the motor to stop increasing in magnitude
before the input power becomes less than the power required, since the area
above Ps would be less than A1.

Figure 5.11 Electric power input to a motor as a function of


torque angle for a suddenly increased load.

Figure 5.12 Electric power input to a motor as a function of torque angle for the
maximum sudden increase of load without loss of stability.
EXAMPLE 5.11 A synchronous motor is receiving 30% of the power which
is capable of receiving from an infinite bus. If the load on the motor is doubled,
calculate the maximum value of d during the swinging of the motor around
its new equilibrium position.
Solution: The load in the motor is 30% and let the initial operating load
angle be d0, i.e. Pe = 0.3Pmax.
244 Electrical Power Systems: Analysis, Security and Deregulation

Figure 5.13

Generally, Pmax sin d0 = Pe = 0.3Pmax


sin d0 = 0.3
d0 = 17.46°
When the load is doubled (Figure 5.13),
Pe = 0.6Pmax
 Pmax sin d1 = 0.6Pmax
sin d1 = 0.6
d1 = 36.87°
To find the maximum value of d (dm)
Using the equal area criterion, A1 = A2
d1

Area A1 = Ú (0.6P
d0
max - Pmax sin d ) d d

= 0.6 Pmax (d1 - d 0 ) + Pmax (cos d )dd10


= 0.6 Pmax (d1 - d 0 ) + Pmax (cos d1 - cos d 0 )
dm

Area A2 = Ú (P
d1
max sin d - 0.6 Pmax ) d d

= - Pmax (cos d m - cos d1 ) - 0.6 Pmax (d m - d1 )


Now A1 = A2
0.6Pmax(d 1 – d 0) + Pmax(cosd 1 – cos d 0)
= –Pmax(cos d m – cos d 1) – 0.6Pmax(d m – d 1)
or 0.6 (d 1 – d 0) + (cos d 1 – cos d 0) = cos d 1 – cos d m – 0.6(d m – d 1)
or –0.6d 0 + cos d 1 – cos d 0 = cos d 1 – cos d m – 0.6d m
or cos d m + 0.6d m = 0.6d 0 + cos d 0
Power System Stability 245

p
cos d m = cos d 0 + 0.6(d 0 - d m ) ¥
180
Ê 0.6 p d m ˆ
= 0.95393 + 0.1828 - Á
Ë 180 ˜¯
cos d m = 1.1368 - 0.01047d m

Solving for dm by trial and error method, we get


d m > d1

d m = 58.15∞

5.7.5 Sudden Three Phase Fault at One End of


the Transmission Line
Let us consider the three-phase short-circuit fault which occurs at the sending
of the transmission line 2 of single machine connected to an infinite system
with double circuit transmission line as shown in Figure 5.14.

Figure 5.14 Three-phase fault at one end of the transmission line.

Prefault condition: Before the occurrence of a fault, both transmission lines


are intact as shown in Figure 5.15. The power angle curve is given by
|E ¢| |V |
Pe1 = sin d = Pmax 1 sin d
XI

Ê X ¥ X TL 2 ˆ
where X I = X d¢ + Á TL1
Ë X TL1 + X TL 2 ˜¯

Figure 5.15 Prefault condition.


246 Electrical Power Systems: Analysis, Security and Deregulation

During fault condition: Upon the occurrence of three phase fault at the
sending of the transmission line 2, the generator gets isolated from the power
system for the purpose of power flow as shown in Figure 5.16. Thus during
the period of fault lasts,
Pe2 = 0

Figure 5.16 During fault condition.

The rotor therefore accelerates and d angle increases. Synchronism will be


lost unless the fault is cleared in time.
Postfault condition: The circuit breakers at the two ends of the faulted line
open at time tc (corresponding to angle dc ), the clearing time, disconnecting
the faulted line as shown in Figure 5.17. The power angle curve is given by
|E ¢| |V |
Pe3 = sin d = Pmax 2 sin d
X III
where XIII = X d + XTL1.

Figure 5.17 Postfault condition.

Obviously, Pmax2 < Pmax1. The rotor now starts to decelerate as shown in
Figure 5.18. The system will be stable if a decelerating area A2 can be found
equal to accelerating area A1 before d reaches the maximum allowable value
dmax. As area A1 depends upon the clearing time tc, the clearing time must be
less than a certain value (critical clearing time) for the system to be stable.
It is to be observed that the equal area criterion helps to determine the critical
clearing angle and not the critical clearing time. Critical clearing time can be
obtained by numerical solution of the swing equation.

5.7.6 Sudden Three Phase Fault at Middle of a


Transmission Line
Consider the fault occurs at middle of a transmission line 2 as shown in
Figure 5.19.
Power System Stability 247

Figure 5.18

Figure 5.19 Fault at middle of transmission line.

Prefault condition: Before the occurrence of a fault, both the lines are
connected as shown in Figure 5.20.

Figure 5.20 Prefault condition.

The power angle curve is given by


|E ¢| |V |
Pe1 = sin d = Pmax 1 sin d
XI
Ê X ¥ X TL 2 ˆ
where, XI = X d¢ + Á TL1
Ë X TL1 + X TL 2 ˜¯

During fault condition: The circuit model of the system during fault is
shown in Figure 5.21. This circuit reduces to that of Figure 5.22 through one
delta to star and one star to delta conversion.
248 Electrical Power Systems: Analysis, Security and Deregulation

Figure 5.21 During fault.

Figure 5.22 Converting reactances from delta to star.

Using delta to star conversion, the circuit becomes as shown in Figure 5.22,
X TL 2
X TL1 ¥
XA = 2 = XB
X TL 2 X TL 2
X TL1 + +
2 2
X TL 2 X TL 2
¥
XF = 2 2
X TL 2 X TL2
X TL1 + +
2 2
Convert star connection to delta connection, the circuit becomes as shown in
Figure 5.23,
( X d¢ + X A ) X F + X F X B + ( X d¢ + X A ) X B
X II =
XF

The power angle curve during fault is given by


|E ¢| |V |
Pe 2 = sin d = Pmax 2 sin d
X II

Figure 5.23 Converting reactances from star to delta.


Power System Stability 249
Postfault condition: In this condition, the faulted transmission line 2 is open
after fault and the circuit is shown in Figure 5.24.
|E ¢| |V |
Pe3 = sin d = Pmax 3 sin d
X III
where XIII = Xd + XTL1.

Figure 5.24 Postfault condition.

Pe1, Pe2, and Pe3 are plotted in Figure 5.25.

Figure 5.25 Power angle curves.

The accelerating area A1 corresponding to d1 is less than area A2, giving


better chance for the stable operation. It is possible to find an area A2 equal
to area A1. As d1 increases, area A1 increases.
To find A1 = A2, d1 increases till d2 = dmax. This case of critical clearing
angle is shown in Figure 5.26.

Figure 5.26 Power angle curves.


250 Electrical Power Systems: Analysis, Security and Deregulation

Applying the equal area criterion, we can write


Area A1 = Area A2
d cr d max

Ú (P m - Pmax 2 sin d ) d d = Ú ( Pmax 3 sin d - Pm ) d d


d0 d cr

Ê P ˆ
where d max = p - sin -1 Á m ˜
Ë Pmax 3 ¯
Integrating, we get
d d
( Pmd + Pmax 2 cos d ) d cr + ( Pmax 3 cos d + Pmd ) d max = 0
0 cr

Pm (d cr - d 0 ) + Pmax 2 (cos d cr - cos d 0 ) + Pm (d max - d cr )


+ Pmax 3 (cos d max - cos d cr ) = 0
Pm (d max - d 0 ) - Pmax 2 cos d 0 + Pmax 3 cos d max
cos d cr =
Pmax 3 - Pmax 2
The critical clearing angle can be calculated from the above equation.
The angles in this equation are in radians. The equation modifies as below if
the angles are in degrees.
p
Pm (d max - d 0 ) - Pmax 2 cos d 0 + Pmax 3 cos d max
cos d cr = 180
Pmax 3 - Pmax 2
EXAMPLE 5.12 Given the circuit as shown in the figure below where a
three-phase fault is applied on one end of a line near circuit breaker CB4.
Find the critical fault clearing angle for clearing the fault with simultaneous
opening of breaker CB2 and CB4. The generator is delivering 1.0 p.u. MW at
the instant preceding the fault. All the p.u. quantities are on the common MVA.

Solution:
Prefault condition Transfer reactance during prefault operation is
0.5 ¥ 0.5
X I = 0.25 + + 0.06 = 0.56
0.5 + 0.5
|E ¢| |V | 1.25 ¥ 1.0
Pe1 = sin d = sin d = 2.232 sin d
XI 0.56
Power System Stability 251
During fault condition The fault occurs at the end of the line 2 or near
bus 2. Therefore during the short circuit fault the circuit separates by the circuit
breaker for finding the transfer reactance as shown in the Figure below. During
the clearing of fault, no power is transferred from the circuit, i.e. Pe2 = 0.

Postfault condition With the opening of the faulted line, say by simultaneous
opening of the circuit breakers CB2 and CB4, the postfault transfer reactance is
X III = 0.25 + 0.5 + 0.06 = 0.81
|E ¢| |V | 1.25 ¥ 1.0
Pe3 = sin d = sin d = 1.543 sin d
X III 0.81
The initial power angle d0 is calculated as
Pm 0 = Pe0 = 1.0 = Pmax 1 sin d 0
Ê 1.0 ˆ
d 0 = sin -1 Á = 26.32∞
        Ë 2.232 ˜¯
and
Ê P ˆ
d max = 180∞ - sin -1 Á m 0 ˜
Ë Pmax 3 ¯
Ê 1.0 ˆ
d max = 180∞ - sin -1 Á = 139.6∞
Ë 1.543 ˜¯
p
Pm (d max - d 0 ) - Pmax 2 cos d 0 + Pmax 3 cos d max
cos d cr = 180
Pmax 3 - Pmax 2
p
1.0(139.6 - 26.62) - 0 + 1.543 cos139.6∞
= 180 = 0.51640
1.543 - 0
d cr = 58.9∞

EXAMPLE 5.13 A 2220 MVA, 24 kV and 60 Hz synchronous machine is


connected to an infinite bus through transformer and double circuit transmission
line, as shown in the following figure. The infinite bus voltage V = 1.0 p.u.
The direct axis transient reactance of the machine is 0.30 p.u., the transformer
reactance is 0.20 p.u., and the reactance of each transmission line is 0.3 p.u., all to
a base of the rating of the synchronous machine. Initially, the machine is delivering
0.8 p.u. real power and reactive power is 0.074 p.u. with a terminal voltage of
1.0 p.u. The inertia constant H = 5 MJ/MVA. All resistances are neglected.
252 Electrical Power Systems: Analysis, Security and Deregulation

(i) A temporary three-phase fault occurs at the sending end of one of


the lines. When the fault is cleared, both lines are intact. Determine
the critical clearing angle and the critical fault clearing time.
(ii) A three-phase fault occurs at the middle of one of the lines, fault is
cleared, and the faulted line is isolated. Determine the critical clearing
angle.

  
Solution:
Given: infinite bus voltage V = 1.0 p.u.
terminal voltage of generator Et =1.0 p.u.
Xd = j0.3
Pm0 = Pe0 = 0.8
The current flowing into the infinite bus is
S* 0.8 - j 0.074
I= = = 0.8 - j 0.074
V *
1.0
The transfer reactance between the internal voltage and the infinite bus before
fault is
0.3 ¥ 0.3
X I = 0.3 + 0.2 + = 0.65
0.3 + 0.3
E ¢ = V + jX I I = 1.0 + j (0.65) (0.8 - j 0.074) = 1.17–26.387∞
(i) Three-phase fault occurs at the sending end of one of the lines
Prefault condition
0.3 ¥ 0.3
X I = 0.3 + 0.2 += 0.65
0.3 + 0.3
|E ¢| |V | 1.17 ¥ 1.0
Pe1 = sin d = sin d = 1.8 sin d
    XI 0.65
During the fault condition: Pe2 = 0
Postfault condition Since both the lines are intact when the fault is
cleared therefore the power angle equation of prefault and postfault
are the same.
|E ¢| |V | 1.17 ¥ 1.0
Pe3 = sin d = sin d = 1.8 sin d
XI 0.65
Power System Stability 253
The initial power angle d0 is calculated as
Pm 0 = Pe0 = 0.8 = Pmax 1 sin d 0
Ê 0.8 ˆ
d 0 = sin -1 Á = 26.388∞ = 0.46055 rad
Ë 1.8 ˜¯
and
Ê P ˆ
d max = 180∞ - sin -1 Á m 0 ˜
Ë Pmax 3 ¯
Ê 0.8 ˆ
d max = 180∞ - sin -1 Á = 153.612∞
          Ë 1.8 ˜¯
p
Pm (d max - d 0 ) - Pmax 2 cos d 0 + Pmax 3 cos d max
cos d cr = 180
Pmax 3 - Pmax 2
p
0.8(153.612 - 26.388) - 0 + 1.8 cos153.612∞
= 180 = 1.48 rad
1.8 - 0
d cr = 84.775∞
The critical clearing time

2 H (d cr - d 0 )
tc =
p fPm
2 ¥ 5(1.48 - 0.46055)
tc = = 0.26 s
p ¥ 60 ¥ 0.8
(ii) Three-phase fault occurs at the middle of one of the lines
Prefault condition From the reactance diagram,
Ê 0.3 ¥ 0.3 ˆ
X T = X I = 0.3 + 0.2 + Á = 0.65
Ë 0.3 + 0.3 ˜¯
|E ¢| |V | 1.17 ¥ 1.0
Pe1 = sin d = sin d = 1.8 sin d
XI 0.65
During the fault condition (fault at the middle of the transmission
line 2)
254 Electrical Power Systems: Analysis, Security and Deregulation

Using delta to star conversion, the circuit becomes

From the reactance diagram,


X T = X II = 1.8
|E ¢| |V | 1.17 ¥ 1.0
Pe 2 = sin d = sin d = 0.65 sin d
X II 1.8
Postfault condition

From the reactance diagram,


X T = X III = 0.8
|E ¢| |V | 1.17 ¥ 1.0
Pe3 = sin d = sin d = 1.4625 sin d
X III 0.8
The initial power angle d0 is calculated as
Pm 0 = Pe0 = 0.8 = Pmax 1 sin d 0
Ê 0.8 ˆ
d 0 = sin -1 Á = 26.388∞ = 0.46055 rad
Ë 1.8 ˜¯

Power System Stability 255
and
Ê P ˆ
d max = 180∞ - sin -1 Á m 0 ˜
Ë Pmax 3 ¯
Ê 0.8 ˆ
d max = 180∞ - sin -1 Á = 146.838∞ = 2.5628 rad
Ë 1.4625 ˜¯
p
Pm (d max - d 0 ) - Pmax 2 cos d 0 + Pmax 3 cos d max
cos d cr = 180
Pmax 3 - Pmax 2
p
0.8(146.838 - 26.388) - 0.65 cos 26.388∞ + 1.4625 cos 146.838∞
= 180
1.4625 - 0.65
= - 0.15356 rad
d cr = 98.834∞

5.8 Solution of Swing Equation


The swing equation, governing the motion of each machine of a system, is
d 2d
= Pa M
dt 2
where d – displacement angle of rotor with respect to a reference axis rotating
at normal speed
M – inertia constant of machine
Pa – accelerating power (difference between the mechanical input and
output after correcting the losses)
t – time.
The solution of swing equation can be done by applying the following
methods.
1. Step by step method.
2. Euler’s method.
3. Modified Euler’s method.
4. Runge–Kutta method.
The step by step method (point by point method) is the most feasible
and widely used way of solving the swing equations. By this method, a good
accuracy can be attained and it is also easy to compute d .
In the step by step method, one or more variables are assumed constant
throughout the short interval of time t, so that as a result of the assumptions
made, the equations can be solved for the changes in the other variables during
the same time interval. Then from the values of the other variables at the end
of the interval, the new values can be calculated for the variables which were
assumed constant. These new values are then used in the next time interval.
256 Electrical Power Systems: Analysis, Security and Deregulation

5.8.1 Step by Step Method-I


It consists of two processes. These are to be carried out alternately.
1. Assume accelerating power Pa to be constant and compute the angular
position and if necessary, the angular velocity (speed) at the end of
the time interval from the knowledge of the position and speeds at the
beginning of the interval.
2. Then from the angular positions and speed, the accelerating power of
each machine is to be calculated. This Pa will be kept constant for the
next time interval and the procedure is to be repeated till the required
final time is reached.

The swing equation is given by


d 2d
M = Pa
dt 2
The above equation is written as,
d 2d Pa
=
dt 2
M
Integrating with respect to t,
dd P
= w = w0 + a t (5.27)
dt M
Integrating once again,
Pa t 2
d = d + w 0t + ◊ (5.28)
M 2
Equations (5.27) and (5.28) respectively provide the speed of the machine (w)
and the angular displacement (d) of the machine with respect to a reference
axis rotating at synchronous speed.
d0 and w0 are the values of d and w at the beginning of the interval.
These equations hold for any instant of time t during the interval in which
Pa is constant.
We are in the need of d and w at the end of the interval.
Let n denote the quantities at the end of the nth interval and n – 1 denote
the quantities at the end of the (n – 1)th interval, which is the beginning of
the nth interval. t is the length of the interval.
Putting t in place of t in Eqs. (5.27) and (5.28) and using the appropriate
subscripts, we obtain the speed and angle at the end of the nth interval as
Dt
w n = w n -1 + Pa (n - 1)
M
(5.29)
Dt
d n = d n -1 + Dt w n -1 + Pa (n - 1)
2M
Power System Stability 257
The increments of speed and angle during the nth interval are
Dt
Dw n = w n - w n -1 = Pa (n - 1)
M
(5.30)
Dt 2
Dd n = d n - d n -1 = Dt w n -1 + Pa (n - 1)
2M
Equation (5.29) or (5.30) is suitable for step by step calculation.
If we are interested only in the angular position and not in the speed,
wn–1 can be eliminated from Eqs. (5.29) and (5.30).
For the preceding interval, we can write
( Dt ) 2
d n -1 = d n - 2 + Dt w n - 2 + Pa ( n - 2) (5.31)
2M
Equation (5.29)–(5.31) gives
Dt 2
d n - d n -1 = d n -1 - d n - 2 + (w n -1 - w n - 2 ) Dt + ( Pa ( n -1) - Pa ( n - 2) )
2M
Putting dn = dn – dn–1; dn–1 = dn–1 – dn–2
wn–1 = wn–1 – wn–2

Dt 2
Dd n = Dd n -1 + Dw n -1 ◊ Dt + ( Pa ( n -1) - Pa ( n - 2) ) (5.32)
2M
But from Eq. (5.30)
Dt
Dw n -1 = Pa ( n - 2)
M
Substituting in Eq. (5.32)
Dt Dt 2
Dd n = Dd n -1 + ◊ Pa ( n - 2) ◊ Dt + ( Pa ( n -1) - Pa ( n - 2) )
M 2M
Dt 2
= Dd n -1 + ( Pa ( n -1) - Pa ( n - 2) )
2M
This equation gives the increment in angle during any interval in terms of the
increment for the previous interval.
Time interval t should be short enough to give the required accuracy.
If it is too short, it will increase the number of calculations to plot a swing
curve, and thus it provides accuracy.

Limitations of step by step method-I


The acceleration during each interval of time t is constant at the value
corresponding to the beginning of the interval.
258 Electrical Power Systems: Analysis, Security and Deregulation

Figure 5.27 shows the true variation of acceleration (a) as a function of


time and the assumed variation of the above method.

Figure 5.27 Variation of acceleration (a) as a function of time.

EXAMPLE 5.14 Consider a 60 Hz machine for which H = 2.7 MJ/MVA


and it is initially operating in steady state with input and output of 1 p.u. and
an angular displacement of 45 electrical degree with respect to an infinite bus
bar. Upon occurrence of a fault, assume that the input remains constant and
the output is given by Pe = d/90°. Calculate and plot swing curve by step by
step method I. Using the time interval t = 0.05 s. Up to t = 1 s.
Solution:
Given Pi = 1 p.u.; S = 1 p.u.
d
Pe =
90∞
d
Also Pa = Pi – Pe = 1 -
90∞
 Pa = 1 – 0.0111d
SH 2.7
M= = = 2.5 ¥ 10 -4
180 f 180 ¥ 60
We have, by method I
Dt
Dw n = Pa ( n -1)
M
Dt 2
Dd n = Dt ◊ w n -1 + Pa ( n -1)
2M
Substituting t and M in the above two equations
0.05 Dt 2
Dw n = Pa ( n -1) ; d n = Dt ◊ w n -1 + Pa ( n -1)
2.5 ¥ 10-4 2M
Dw n = 200 Pa ( n -1) Dd n = 0.05 w n -1 + 5 Pa ( n -1)

Power System Stability 259

ts Pe Pa w w 0.05 w 5Pa d deg. d deg.


0+ 0.5 0.5 100 0 0 2.5 2.5 45
0.05 0.528 0.472 94.55 100 5 2.36 7.36 47.5
0.1 0.6089 0.391 78.2 194.55 9.7275 1.955 11.6825 54.86
0.15 0.7386 0.2614 52.2757 272.75 13.6375 1.307 14.9445 66.5425
0.2 0.9045 0.0955 19.1 325.0257 16.2513 0.4775 16.7288 81.487
0.25 1.0902 –0.0902 –18.04 344.13 17.21 –0.45 16.76 98.2158
0.3 1.28 –0.28 –55.24 326.09 16.3 –1.4 14.9 114.97
0.35 1.441 –0.44 –88.31 270.85 13.54 –2.2 11.34 129.87
0.4 1.57 –0.57 –113.49 182.54 9.13 –2.85 6.28 141.21
0.45 1.64 –0.64 –127.43 69.05 3.45 –3.2 0.25 147.49
0.5 1.64 –0.64 –127.98 –58.38 –2.92 –3.2 –6.12 147.74
0.55 1.57 –0.57 –114.4 –186.36 –9.32 –2.85 –12.17 141.62
0.6 1.44 –0.44 –88 –300.76 –15.04 –2.2 –17.24 129.45
0.65 1.25 –0.25 –49.11 –388.76 –19.44 –1.25 –20.69 112.21
0.7 1.02 –0.02 –3.17 –437.87 –21.89 –0.1 –21.99 91.52
0.75 0.77 0.23 46 –441.04 –22.05 1.15 –20.9 69.53
0.8 0.54 0.46 92 –395.04 –19.75 2.3 –17.45 48.63
0.85 0.35 0.65 130.75 –303.04 –15.15 3.25 –11.9 31.18
0.9 0.21 0.79 157.2 –172.26 –8.61 3.93 –4.68 19.28
0.95 0.16 0.84 167.59 –15.06 –0.75 4.2 3.45 14.6
1 0.2 0.8 159.94 152.53 7.63 4 11.63 18.05
29.68

The variation of d with time t is plotted as shown in the figure below


From the swing curve, we see that the value of d increases and then
decreases, therefore the system is stable. A more accurate approach would be to
assume the average value of acceleration over the interval t, by using the value
of acceleration at the middle of the interval. It is adopted in step by step method II.
260 Electrical Power Systems: Analysis, Security and Deregulation

MATLAB program—solution of swing equation using step by


step method – I
Del_delta=0; i=1;
f=input(‘Enter the frequency :’);
H=input(‘Enter the value of inertia constant :’);
delta=input(‘Enter initial displacement angle :’);
Pi=input(‘Enter initial steady state power :’);
M=H/(180*f); y=0.05/M;yy=0.05^2/(2*M);Pu=delta/90;
del_int=delta*pi/180; Pa=1-Pu;
disp(‘Time Pe Pa 5*Pa Del_ang Del’)
for t=0:.05:1.05
if t==0
w=0;
del_w=y*Pa;
del_del=0.05*w+yy*Pa;
fprintf(‘%g-’,t);
disp([Pi Pa yy*Pa del_del delta]);
fprintf(‘%g+’,t);
disp([Pu Pa yy*Pa del_del delta])
else
w=w+del_w;
delta=delta+del_del;
Pu=delta/90;
Pa=1-Pu;
del_w=y*Pa;
del_del=0.05*w+yy*Pa;
fprintf(‘%g’,t);
disp([Pu Pa yy*Pa del_del delta])
end
time(i)=t;
del(i)=delta;
delta=delta+Del_delta;
i=i+1;
end
plot(time,del);
title(‘SWING CURVE’);
xlabel(‘t, sec’);
ylabel(‘delta, elec. deg’);
Results:
Enter the value of inertia constant: 2.7
Enter the initial displacement angle: 45
Enter the initial steady state power: 1
Power System Stability 261

Time Pe Pa 5*Pa Del_ang Del


0+ 0.5 0.5 2.5 2.5 45
0.05 0.52778 0.47222 2.3611 7.3611 47.5
0.1 0.60957 0.39043 1.9522 11.674 54.861
0.15 0.73928 0.26072 1.3036 14.93 66.535
0.2 0.90517 0.09483 0.47413 16.708 81.466
0.25 1.0908 –0.0908 –0.4541 16.728 98.173
0.3 1.2767 –0.2767 –1.3834 14.89 114.9
0.35 1.4421 –0.4421 –2.2107 11.296 129.79
0.4 1.5676 –0.5677 –2.8382 6.2475 141.09
0.45 1.6371 –0.6371 –3.1853 0.22392 147.34
0.5 1.6395 –0.6396 –3.1977 –6.1591 147.56
0.55 1.5711 –0.5711 –2.8556 –12.212 141.4
0.6 1.4354 –0.4354 –2.1771 –17.245 129.19
0.65 1.2438 –0.2438 –1.219 –20.641 111.94
0.7 1.0145 –0.0145 –0.0723 –21.933 91.301
0.75 0.77076 0.22924 1.1462 –20.859 69.369
0.8 0.539 0.461 2.305 –17.408 48.51
0.85 0.34558 0.65442 3.2721 –11.831 31.102
0.9 0.21413 0.78587 3.9293 –4.6291 19.272
0.95 0.1627 0.8373 4.1865 3.4868 14.643
1 0.20144 0.79856 3.9928 11.666 18.13
262 Electrical Power Systems: Analysis, Security and Deregulation

5.8.2 Step by Step Method-II


Acceleration is assumed constant from the middle of one interval to the middle
of the next interval.
For example, consider Figure 5.28(a) where an–1 is the equivalent constant
value of acceleration from t = (n – (3/2))t to t = (n – (1/2))t.

Figure 5.28 Step by step method-II—variation of a, w and d.


Power System Stability 263
an is constant from t = (n – (1/2))t to t = (n + (1/2))t and so on. In the
region of constant acceleration an–1, there is an increment of angular velocity
from wn–3/2 to wn–1/2, i.e. [from Figure 5.28(b)]
wn–1/2 = wn–3/2 + wn–1/2 (5.33)
where wn–1/2 is the increment in angular velocity over a time interval t,
due to acceleration an–1.
wn–1/2 = an–1 · t (5.34)

d 2d n -1 Pa ( n -1)
As an–1 = 2
=
dt M (5.35)
Pa ( n -1)
Thus, wn–1/2 = ◊ Dt
M
Again, the angular velocity, wn–1/2 remains constant for t = (n – 1)t, during
the nth interval. From Figure 5.28(c), the displacement angle dn–1 increases to
dn over this interval by an amount dn.

dn = dn–1 + dn (5.36)


dn = wn–1/2 · t (5.37)
Using Eq. (5.33)
wn–1/2 = wn–3/2 + wn–1/2
Hence
dn = (wn–3/2 + wn–1/2)t
(5.38)
= wn–3/2 t + an–1 t2   by using Eq. (5.34)

From Eq. (5.37)


dn–1 = (wn–3/2)t (5.39)

Substituting Eq. (5.39) in Eq. (5.38), we get

dn = dn–1 + an–1 t 2

Substituting an–1 = Pa(n–1) /M, we get


Pa ( n -1)
Dd n = Dd n -1 + Dt 2 (5.40)
M
In the method II, it is not necessary to calculate w, unless it is specifically
required.
Using Eq. (5.40), Pa is to be found step by step over a number of intervals
of time, from which the increment in displacement angle can be calculated as
the inertia constant is known.
264 Electrical Power Systems: Analysis, Security and Deregulation

To begin with, the power output at t = 0– and 0+ are averaged out and
the average value of Pa is determined. Then d is calculated from the value
during the preceding interval.
EXAMPLE 5.15 Consider a 60 Hz machine for which H = 2.7 MJ/MVA
and it is initially operating in steady state with input and output of 1 p.u.
and an angular displacement of 45 electrical degree with respect to an
infinite bus bar. Upon occurrence of a fault, assume that the input remains
constant and the output is given by Pe = d/90°. Calculate and plot the swing
curve by the step-by-step method II. Using the time interval t = 0.05 s.
Up to t = 1 s. Step-by-step method-II using the time interval t = 0.05 s and
upto t = 1 s.
Solution:
H = 2.7 MJ/MVA
At time t = 0–,    Pm = Pe = 1 p.u.
At time t = 0+, Pm = 1 p.u., Pe = d/90°
SH
M = = 2.5 ¥ 10-4 p.u.
180 f
d
Pa = 1 -
90∞
2
Dt 0.052
= = 10
M 2.5 ¥ 10 -4
Now the equation for dn is given by

Dt 2
Dd n = Dd n -1 + Pa ( n -1)
M
= Dd n -1 + 10 Pa ( n -1)

t (s) Pe Pa 10Pa d deg. d deg.



0 1 0 — — —
0+ 0.5 0.5 — — 45°
0av — 0.25 2.5 — —
2.5
0.05 0.5277 0.472 4.72 — 47.5
7.22
0.1 0.608 0.392 3.92 — 54.72
11.14
0.15 0.732 0.2682 2.682 — 65.86
13.822
(Contd...)
Power System Stability 265
(Contd...)

t (s) Pe Pa 10Pa d deg. d deg.


0.2 0.8854 0.1146 1.1464 — 79.68
14.9764
0.25 1.05174 –0.0517 –0.5173 — 94.6564
14.459
0.3 1.212 –0.2124 –2.1239 — 109.115
12.335
0.35 1.3495 –0.3495 –3.495 — 121.45
8.84
0.4 1.448 –0.448 –4.48 — 130.29
4.36
0.45 1.496 –0.496 –4.96 — 134.65
–0.6
0.5 1.489 –0.489 –4.89 — 134.05
–5.49
0.55 1.428 –0.428 –4.28 — 128.56
–9.77
0.6 1.32 –0.32 –3.2 — 118.79
–12.97
0.65 1.176 –1.1758 –1.757 — 105.82
–14.727
0.7 1.012 –0.012 –0.12 — 91.093
–14.877
0.75 0.847 0.157 1.53 — 76.216
–13.347
0.8 0.699 0.301 3.01 — 62.869
–10.337
0.85 0.584 0.416 4.16 — 52.532
–6.177
0.9 0.515 0.485 4.85 — 46.355
–1.327
0.95 0.5 0.5 5 45.028
3.673
1 48.701
266 Electrical Power Systems: Analysis, Security and Deregulation

From the swing curve, we see that the value of d increases and then
decreases, which makes the system stable.
Note: In the method-II, the average of the accelerating powers was calculated
from the prefault (0–) and during the fault (0+) powers, and used in the further
iterations. Similarly, if the fault is cleared by some means at any intermediate
time t, the same method of calculation of average power is to be adopted for
during the fault (t –) and the postfault (t+) time periods, and these values are
to be used in further calculations.

MATLAB program—Solution of swing equation using step by


step method–II
Del_delta=0; i=1;
f=input(‘Enter the frequency :’);
H=input(‘Enter the value of inertia constant :’);
delta=input(‘Enter initial displacement angle :’);
Pi=input(‘Enter initial steady state power :’);
M=H/(180*f); y=0.05/M;yy=0.05^2/(M);Pu=delta/90;
del_int=delta*pi/180; Pa=1-Pu;
disp(‘Time Pe Pa 10*Pa Del_ang Del’)
for t=0:.05:1.05
if t==0
w=0;
del_w=y*Pa;
del_del=0.05*w+yy*Pa;
fprintf(‘%g-’,t);
disp([Pi Pa del_w del_del delta]);
fprintf(‘%g+’,t);
Power System Stability 267
disp([Pu Pa del_w del_del delta])
Pa=Pa/2;
del_w=y*Pa;
del_del=0.05*w+yy*Pa;
fprintf(‘%gavg’,t);
disp([Pu Pa yy*Pa del_del delta])
else
w=w+del_w;
delta=delta+del_del;
Pu=delta/90;
Pa=1-Pu;
del_w=y*Pa;
del_del=0.05*w+yy*Pa;
fprintf(‘%g’,t);
disp([Pu Pa yy*Pa del_del delta])
end
time(i)=t;
del(i)=delta;
delta=delta+Del_delta;
i=i+1;
end
plot(time,del);
title(‘SWING CURVE’);
xlabel(‘t, sec’);
ylabel(‘delta, elec. deg’);
Results:
Enter the value of inertia constant : 2.7
Enter the initial displacement angle : 45
Enter the initial steady state power : 1
Time Pe Pa 10*Pa Del_ang Del
0– 1 0.5 — — 45
0+ 0.5 0.5 — — 45
0avg — 0.25 2.5 — 45
0.05 0.52778 0.47222 4.7222 7.2222 47.5
0.1 0.60802 0.39198 3.9198 11.142 54.722
0.15 0.73182 0.26818 2.6818 13.824 65.864
0.2 0.88542 0.11458 1.1458 14.97 79.688
0.25 1.0517 –0.0517 –0.5174 14.452 94.657
0.3 1.2123 –0.2123 –2.1233 12.329 109.11
0.35 1.3493 –0.3493 –3.4931 8.8356 121.44
0.4 1.4475 –0.4474 –4.4749 4.3607 130.27
(Contd...)
268 Electrical Power Systems: Analysis, Security and Deregulation

(Contd...)

Time Pe Pa 10*Pa Del_ang Del


0.45 1.4959 –0.4959 –4.9594 –0.5986 134.63
0.5 1.4893 –0.4892 –4.8929 –5.4915 134.04
0.55 1.4283 –0.4282 –4.2827 –9.7742 128.54
0.6 1.3197 –0.3196 –3.1967 –12.971 118.77
0.65 1.1755 –0.1755 –1.7555 –14.726 105.8
0.7 1.0119 –0.0119 –0.1192 –14.846 91.073
0.75 0.84697 0.15303 1.5303 –13.315 76.227
0.8 0.69902 0.30098 3.0098 –10.306 62.912
0.85 0.58452 0.41548 4.1548 –6.1507 52.606
0.9 0.51618 0.48382 4.8382 –1.3124 46.456
0.95 0.50159 0.49841 4.9841 3.6717 45.143
1 0.54239 0.45761 4.5761 8.2478 48.815

From the swing curve, the value of d increases and then decreases, which
makes the system stable.
EXAMPLE 5.16 A 20 MVA, 3-f, 50 Hz generator delivers rated power at
unity power factor via a double circuit transmission line to an infinite bus
bar. The generator unit has a kinetic energy of 2.5 MJ/MVA at rated speed.
Its X d¢ = 0.3 p.u. The transformer circuits have negligible resistances and
each has a reactance of 0.3 p.u. on a 20 MVA base. The voltage behind the
transient reactance is 1.05 p.u. and the voltage of the infinite bus is 1 p.u.
A 3-f short circuit occurs at the middle of one of the transformer circuits.
It involves ground (a) What is the initial displacement angle of the machine?
The fault is cleared in 0.4 s by simultaneous opening of CBs at both ends of
Power System Stability 269
the faulted transmission line. (b) Calculate and plot the swing curve for the
system and ascertain whether the system is stable or not.
Take t = 0.05 s and tmax = 1 s.
Solution:
SH 1 ¥ 2.5
M = = = 2.778 ¥ 10 -4
180 f 180 ¥ 50
MW 20 ¥ 1 20
Pi = = = = 1 p.u.
MVA 20 20
The power angle equation for the three conditions should be taken into account.
1. Prefault;   2. During fault;   3. Postfault

Prefault condition
0.3 ¥ 0.3
X1 = X AB = 0.3 + = 0.3 + 0.15
0.3 + 0.3
X1 = 0.45 p.u.
During fault
X2 = XAB is obtained by converting the star connection to delta connection.

0.3 ¥ 0.3
\ X 2 = X AB = 0.3 + 0.3 + = 0.6 + 0.6
0.15
X 2 = 1.2 p.u.
Postfault The faulty transmission line is made out of service by simultaneous
opening of the circuit breakers at both ends and hence the faulty line is removed.
270 Electrical Power Systems: Analysis, Security and Deregulation

X3 = XAB = 0.3 + 0.3


= 0.6 p.u.

1.05 ¥ 1
Prefault power, Pmax 1 = = 2.333 p.u.
0.45
1.05 ¥ 1
During fault power, Pmax 2 = = 0.875 p.u.
1.2
1.05 ¥ 1
Postfault power, Pmax 3 = = 1.75 p.u.
0.6
We know that
Dt 2
dn = Dd n -1 + Pa ( n -1)
M
Dt 2 (0.05) 2
 = =9
M 2.7778 ¥ 10-4
 dn = dn–1 + 9Pa(n–1)
Pi = 1 p.u.
Pa = 1 – Pe
In general, Pe = Pmax sind
For the prefault condition

Pmax 1 sin d = Pm
2.333 sin d = 1
1
sin d = = 0.4286
2.333
 d = 25.3808
For the during fault condition, Pe = Pmax2 sin d
For the postfault condition, Pe = Pmax3 sin d
Power System Stability 271
The calculation of d versus t is done in the following table.
t (s) Pmax sind Pe = Pa 9Pa d d degree
Pmax sin d degree
0–­ 2.333 0.4286 1 — — — 25.3808
0+ 0.875 0.4286 0.375 0.625 — — 25.3808
0av — — — 0.3125 2.8125 — —
2.8125
0.05 0.875 0.4724 0.4134 0.5866 5.2795 — 28.1933
8.0920
0.1 0.875 0.5918 0.5178 0.4822 4.3395 — 36.2853
12.4315
0.15 0.875 0.7515 0.6575 0.3425 3.0823 — 48.7168
15.5138
0.2 0.875 0.9006 0.7880 0.2120 1.9082 — 64.2306
17.4220
0.25 0.875 0.9894 0.8657 0.1343 1.2084 — 81.6526
18.6304
0.3 0.875 0.9839 0.8609 0.1391 1.2515 — 100.2830
19.8819
0.35 0.875 0.8646 0.7565 0.2435 2.1914 — 120.1649
22.0733
0.4– 0.875 0.6124 0.5358 — — — 142.2382
0.4+ 1.75 0.6124 1.0717 — — — 142.2382
0.4av — — 0.8038 0.1962 1.7658 —
23.8391
0.45 1.75 0.2406 0.4211 0.5789 5.2104 166.0773
29.0495
0.5 1.75 –0.2610 –0.4567 1.4567 13.1100 195.1268
42.1595
0.55 1.75 –0.8414 –1.4724 2.4724 22.2518 237.28
64.4113
0.6 1.75 –0.8508 –1.489 2.489 22.4006 301.6976
86.8119
0.65 1.75 0.4773 0.8353 0.1647 1.4825 388.5095
88.2944
0.7 1.75 0.8926 1.5620 –0.5620 –5.0577 476.8039
83.2367
0.75 1.75 –0.3427 –0.5997 1.5997 14.3973 560.0406
97.634
0.8 1.75 –0.8856 –1.5498 2.5498 22.9482 657.67
120.5822
(Contd...)
272 Electrical Power Systems: Analysis, Security and Deregulation

(Contd...)

t (s) Pmax sind Pe = Pa 9Pa d d degree


Pmax sin d degree
0.85 1.75 0.8504 1.4882 –0.4882 –4.394 778.2568
116.1882
0.9 1.75 0.0968 0.1694 0.8306 7.4754 894.445
123.6636
0.95 1.75 –0.8821 –1.5436 2.5436 22.8924 1018.1086
146.556
1.0 1.75 0.9957 1.7424 –0.7424 –6.6818 1164.6646
139.8742
1.05 1304.5388

The variation of d with time t is plotted as shown below.

(a) Initial displacement angle d0 = 25.3808.


(b) From the swing curve, the displacement angle d goes on increasing.
Hence the system is unstable.

MATLAB program—solution of swing equation by considering


critical clearing time using step by step method–II
clc;
clear;format short g;
Pm=input(‘Enter the power limits for prefault,during
fault,postfault conditions :’);
Del_delta=0; i=1;
Power System Stability 273
f=input(‘Enter the frequency :’);
tc=input(‘Enter the Fault Clearing Time :’);
H=input(‘Enter the value of inertia constant :’);
delta=input(‘Enter initial displacement angle :’);
Pi=input(‘Enter initial steady state power :’);
M=H/(180*f); y=0.05^2/M;
k1=Pm(2)/Pm(1);k2=Pm(3)/Pm(1);
del_int=delta*pi/180;
del_max=pi-asin(sin(del_int)/k2);
del_cri=acos(((del_max-del_int)*sin(del_int)-k1*cos(del_
int)+k2*cos(del_max))/(k2-k1));
cri_t=(2*M*(del_cri-del_int)/Pi)^(1/2);
del_cri=(del_cri*180)/pi;
disp(‘Time P_max Sin(del) Pe Pa 9*Pa Del_ang Delta\n’)
for t=0:.05:1.05
delta=delta*pi/180;
if t==0
Paminus=Pi-1;
Paplus=Pi-Pm(2)*sin(delta);
Paavg=(Paminus+Paplus)/2;Pma=(Pm(1)+Pm(2))/2;sd=sin(delta);
Pa=Paavg;delta=(delta*180)/pi;
fprintf(‘%g-’,t);
disp([Pm(1) sd Pi Paminus (y*Paminus) Del_delta delta]);
fprintf(‘%g+’,t);
disp([Pm(2) sd Pm(2)*sd Paplus (y*Paplus) Del_delta delta])
fprintf(‘%gavg’,t);
disp([Pma sd Pma*sd Pa (y*Pa) Del_delta delta])
end
if t==tc
Paminus=Pi-Pm(2)*sin(delta);
Paplus=Pi-Pm(3)*sin(delta);
Paavg=(Paminus+Paplus)/2;Pma=(Pm(3)+Pm(2))/2;sd=sin(delta);
Pa=Paavg;delta=(delta*180)/pi;
fprintf(‘%g-’,tc);
disp([Pm(2) sd Pm(2)*sd Paminus (y*Paminus) Del_delta delta])
fprintf(‘%g+’,tc);
disp([Pm(3) sd Pm(3)*sd Paplus (y*Paplus) Del_delta delta])
fprintf(‘%gavg’,tc);
disp([Pma sd Pma*sd Pa (y*Pa) Del_delta delta])
end
if t>0 && t<tc
Pa=Pi-Pm(2)*sin(delta);sd=sin(delta);delta=(delta*180)/pi;
disp([t Pm(2) sd Pm(2)*sd Pa (y*Pa) Del_delta delta])
end
274 Electrical Power Systems: Analysis, Security and Deregulation

if (t>tc)
Pa=Pi-Pm(3)*sin(delta);sd=sin(delta);delta=(delta*180)/pi;
disp([t Pm(3) sd Pm(3)*sd Pa (y*Pa) Del_delta delta])
end
Del_delta=Del_delta+(y*Pa);
time(i)=t;
del(i)=delta;
delta=delta+Del_delta;
i=i+1;
end
critical_clearing_angle=del_cri
critical_clearing_time=cri_t
plot(time,del);
title(‘SWING CURVE’);
xlabel(‘t, sec’);
ylabel(‘delta, elec. deg’);
Results:
Enter the power limits for prefault, during fault, and postfault conditions:
[2.333 .875 1.75]
Enter the frequency: 50
Enter the fault clearing time: 0.4
Enter the value of inertia constant: 2.5
Enter the initial displacement angle: 25.3808
Enter the initial steady state power: 1.

Time P_max sin(del) Pe Pa 9*Pa Del_ang Delta


0– 2.333 0.42863 1 — — — 25.381
0+ 0.875 0.42863 0.37505 0.62495 — — 25.381
0avg — — — 0.31247 2.8123 — —

0.05 0.875 0.47244 0.41339 0.58661 5.2795 2.8123 28.193


0.1 0.875 0.5918 0.51782 0.48218 4.3396 8.0918 36.285
0.15 0.875 0.75145 0.65752 0.34248 3.0823 12.431 48.716
0.2 0.875 0.90055 0.78798 0.21202 1.9082 15.514 64.23
0.25 0.875 0.9894 0.86573 0.13427 1.2084 17.422 81.652
0.3 0.875 0.98394 0.86095 0.13905 1.2515 18.63 100.28
0.35 0.875 0.86459 0.75652 0.24348 2.1913 19.882 120.16
0.4– 0.875 0.6124 0.53585 0.46415 4.1774 22.073 142.24
0.4+ 1.75 0.6124 1.0717 –0.0717 –0.6452 22.073 142.24
0.4avg 1.3125 0.6124 0.80377 0.19623 1.7661 22.073 142.24
0.45 1.75 0.24063 0.42111 0.57889 5.21 23.839 166.08
0.5 1.75 –0.26093 –0.4566 1.4566 13.11 29.049 195.13
(Contd...)
Power System Stability 275
(Contd...)
Time P_max sin(del) Pe Pa 9*Pa Del_ang Delta
0.55 1.75 –0.84136 –1.4724 2.4724 22.251 42.159 237.28
0.6 1.75 –0.85086 –1.489 2.489 22.401 64.41 301.69
0.65 1.75 0.47725 0.83519 0.16481 1.4833 86.811 388.51
0.7 1.75 0.89258 1.562 –0.5620 –5.0582 88.295 476.8
0.75 1.75 –0.34263 –0.5996 1.5996 14.396 83.237 560.04
0.8 1.75 –0.88564 –1.5499 2.5499 22.949 97.633 657.67
0.85 1.75 0.85037 1.4881 –0.4881 –4.3933 120.58 778.25
0.9 1.75 0.096883 0.16954 0.83046 7.4741 116.19 894.44
0.95 1.75 –0.8821 –1.5437 2.5437 22.893 123.66 1018.1
1 1.75 0.99566 1.7424 –0.7424 –6.6816 146.56 1164.7
1.05 1.75 –0.70131 –1.2273 2.2273 20.046 139.87 1304.5

critical clearing angle = 98.963


critical clearing time = 0.026711

Initial displacement angle, d0 = 25.3808.


From the swing curve, the displacement angle d goes on increasing.
Hence the system is unstable.
EXAMPLE 5.17 A 25 MVA, 60 Hz, water wheel generator delivers
20 MW over a double circuit transmission line to a large metropolitan system
which may be regarded as an infinite bus. The generating unit has a kinetic
energy of 2.76 MJ/MVA at rated speed. The direct axis transient reactance of
276 Electrical Power Systems: Analysis, Security and Deregulation

the generator is 0.3 p.u. The transmission circuits have negligible resistances
and each has a reactance of 0.2 p.u. on a 25 MVA base. The voltage behind
the transient reactance of the generator is 1.03 p.u. and the voltage of the
metropolitan system is 1.0 p.u. A three-phase short circuit occurs at the middle
of the transmission line circuit and is cleared in 0.4 s by the simultaneous
opening of the circuit breaker at both ends of the line. Calculate and plot the
swing curve of the generator for 1 s.
Solution:

The net reactance for the prefault condition is given by


0.2 ¥ 0.2
X1 = X AB = 0.3 + = 0.4 p.u.
0.2 + 0.2
During the fault condition

0.3 ¥ 0.2 + 0.2 ¥ 0.1 + 0.1 ¥ 0.3


X 2 = X AB =
0.1
0.06 + 0.02 + 0.03
= = 1.1 p.u..
0.1
X 2 = 1.1 p.u.

During the postfault condition


X3 = XAB = 0.3 + 0.2 = 0.5 p.u.
Power System Stability 277

1.03 ¥ 1
Prefault power, Pmax 1 = = 2.58 p.u.
0.4
1.03 ¥ 1
During fault power, Pmax 2 = = 0.936 p.u.
1.1
1.03 ¥ 1
Postfault power, Pmax 3 = = 2.06 p.u.
0.5
SH 1 ¥ 2.76
M = = = 2.56 ¥ 10 -4
180 f 180 ¥ 60

20
Output, Pe = = 0.8 p.u.
25
In general Pe = Pmaxsind
For the prefault condition,
Pe = Pmax 1 sin d
0.8 = 2.58 sin d
0.8
sin d =
2.58
0.8
d = sin -1 = 18.2∞
2.58
For the during fault condition, Pe = Pmax2 sin d
For the postfault condition, Pe = Pmax3 sin d
Immediately after the fault, d remains unchanged momentarily but output
changes.
Pe = Pmax 2 sin d = 0.936 sin 18.2∞
= 0.292 p.u.
Dt 2
Dd n = Dd n -1 + Pa ( n -1)
M
(0.05) 2
= Dd n -1 + Pa ( n -1)
2.564 ¥ 10-4
= Dd n -1 + 9.76 Pa ( n -1)

278 Electrical Power Systems: Analysis, Security and Deregulation

The calculation of d versus t is done in the following table.


Pa = Pi – Pe = 0.8 – Pe
t (s) Pmax sind Pu Pa 9.76Pa d d degree
degree
0– 2.58 0.31 0.8 0 — — 18.1
0+ 0.936 0.31 0.29 0.51 — — 18.1
0av — — — 0.255 2.5 — —
2.5
0.05 0.936 0.352 0.33 0.47 4.6 — 20.6
7.1
0.1 0.936 0.465 0.435 0.365 3.562 — 27.7
10.662
0.15 0.936 0.621 0.581 0.219 2.138 — 38.362
12.800
0.2 0.936 0.779 0.729 0.071 0.692 — 51.162
13.492
0.25 0.936 0.904 0.846 –0.046 –0.448 — 64.654
13.044
0.3 0.936 0.977 0.915 –0.115 –1.118 — 77.698
11.926
0.35 0.936 1.0 0.936 –0.136 –1.327 — 89.624
10.599
0.4– 0.936 0.984 0.921 –0.121 — — 100.223
0.4+ 2.06 0.983 2.025 –1.125 — — 100.223
0.4av — — — –0.623 –6.568 — —
4.031
0.45 2.06 0.969 1.997 –1.197 –11.679 — 104.254
–7.648
0.5 2.06 0.993 2.046 –1.246 –12.164 — 96.606
–19.812
0.55 2.06 0.974 2.006 –1.206 –11.766 — 76.794
–31.578
0.6 2.06 0.710 1.462 –0.662 –6.462 — 45.216
–38.040
0.65 2.06 0.125 0.257 0.543 5.297 — 7.176
–32.743
0.7 2.06 –0.432 –0.890 1.690 16.494 — –25.567
–16.249
0.75 2.06 –0.667 –1.373 2.173 21.213 — –41.816
4.964
(Contd...)
Power System Stability 279
(Contd...)
t (s) Pmax sind Pu Pa 9.76Pa d d degree
degree
0.8 2.06 –0.600 –1.235 2.035 19.866 — –36.852
24.83
0.85 2.06 –0.208 –0.429 1.229 11.996 –12.022
36.826
0.9 2.06 0.420 0.864 –0.064 –0.627 — 24.804
36.199
0.95 2.06 0.875 1.802 –1.002 –9.777 — 61.003
26.422
1.0 2.06 0.999 2.058 –1.258 –12.277 — 87.425

The variation of d with time t is plotted as shown below.

If we plot the swing curve, the value of d increases and then decreases,
this makes the system stable.
Results:
Enter the power limits for prefault, during fault, and postfault conditions:
[2.58 .936 2.06]
Enter the frequency: 60
Enter the fault clearing time: 0.4
Enter the value of inertia constant: 2.76
Enter the initial displacement angle: 18.2
Enter the initial steady state power: 0.8.
280 Electrical Power Systems: Analysis, Security and Deregulation

t (s) Pmax sind Pe Pa 9*Pa d degree d degree


0– 2.58 0.31233 0.8 0 — 0 18.2
0+ 0.936 0.31233 0.29235 0.50765 — 0 18.2
0avg — — — 0.255 2.5 0 18.2
0.05 0.936 0.33717 0.3156 0.4844 4.7387 1.5048 19.705
0.1 0.936 0.43756 0.40956 0.39044 3.8195 6.2436 25.948
0.15 0.936 0.58795 0.55032 0.24968 2.4425 10.063 36.012
0.2 0.936 0.74915 0.70121 0.098792 0.96644 12.506 48.517
0.25 0.936 0.88286 0.82636 –0.02636 –0.2578 13.472 61.989
0.3 0.936 0.96684 0.90496 –0.10496 –1.0268 13.214 75.204
0.35 0.936 0.99896 0.93503 –0.13503 –1.3209 12.187 87.391
0.4– 0.936 0.98963 0.9263 –0.1263 –1.2355 10.867 98.257
0.4+ 2.06 0.98963 2.0386 –1.2386 –12.117 10.867 98.257
0.4avg 1.498 0.98963 1.4825 –0.68247 –6.6763 10.867 98.257
0.45 2.06 0.97649 2.0116 –1.2116 –11.852 4.1902 102.45
0.5 2.06 0.99651 2.0528 –1.2528 –12.256 –7.6622 94.785
0.55 2.06 0.96532 1.9886 –1.1886 –11.627 –19.918 74.867
0.6 2.06 0.6861 1.4134 –0.61336 –6.0003 –31.545 43.322
0.65 2.06 0.10065 0.20733 0.59267 5.7978 –37.546 5.7764
0.7 2.06 –0.4379 –0.9021 1.7021 16.651 –31.748 –25.971
0.75 2.06 –0.6569 –1.3533 2.1533 21.065 –15.097 –41.068
0.8 2.06 –0.575 –1.1845 1.9845 19.414 5.9685 –35.099
0.85 2.06 –0.1687 –0.3477 1.1477 11.228 25.382 –9.7174
0.9 2.06 0.45231 0.93176 –0.13176 –1.289 36.61 26.892
0.95 2.06 0.88468 1.8225 –1.0225 –10.002 35.321 62.213
1 2.06 0.99907 2.0581 –1.2581 –12.307 25.318 87.531
1.05 2.06 0.98312 2.0252 –1.2252 –11.986 13.011 100.54
critical clearing angle = 137.85
critical clearing time = 0.036526

In the swing curve, the value of d increases and then decreases, this makes
the system stable.
Power System Stability 281
5.8.3 Euler’s Method
The Euler’s method is the simplest and the least accurate of all numerical
methods. It is presented here because of its simplicity. By studying this method,
we will be able to grasp the basic ideas involved in the numerical solutions of
one-dimensional equation (ODE) and can easily understand the comparatively
complex method such as the Runge–Kutta procedure.
Let us consider the first order differential equation
dx
= f ( x, t ) (5.41)
dt
Figure 5.29 illustrate the principles of applying the Euler’s method at initial
condition x = x0 at t = t0.

Figure 5.29 Graphical interpretation of Euler’s method.

At x = x0, t = t0, we can approximate the curve representing the true


solution by its tangent having a slope
dx
= f ( x0 , t0 ) (5.42)
dt x = x0

For a small increment in t denoted by t, the increment in x is given by x


Therefore,
dx
Dx = ◊ Dt (5.43)
dt x = x0

dx
where is the slope of the curve at (t0, x0), which can be determined
dt x = x0
from Eq. (5.42).Thus, the value of x at t = t1 = t0 + t is given by

dx
x1 = x0 + Dx = x0 + ◊ Dt (5.44)
dt x = x0
282 Electrical Power Systems: Analysis, Security and Deregulation

The Euler’s method is equivalent to using the first two terms of the Taylor
series expansion for x around the point (t0, x0)

Dt 2 Dt 3
x1 = x0 + Dt ( x0 ) + ( 
x0 ) + (
x0 ) +  (5.45)
2! 3!
After using the Euler’s technique for determining x = x1 corresponding to
t = t1, we can take another short time step t and determine x2 corresponding
to t2 = t1 + t as follows:
dx
x2 = x1 + ◊ Dt (5.46)
dt x = x1
The subsequent values of x can be similarly determined. Hence, the computational
algorithm is
dx
xi +1 = xi + ◊ Dt (5.47)
dt x = xi

By applying the above algorithm successively, the values of x can be determined


corresponding to different values of t.
The method considers only the first derivative of x and is, therefore,
referred to as a first order method. For sufficient accuracy at each step t has
to be small. This will increase round-off errors, and the computational effort
required will be very high.

5.8.4 Modified Euler’s Method


The standard Euler’s method results in inaccuracies because it uses the derivative
at the beginning of the interval though it applied throughout the interval.
The slope is constant over the entire interval t causing the points to fall
below the curve. The above problem is solved by the modified Euler’s method.
This method can be obtained by calculating the slope both at the beginning
and the end of the interval, and then averaging these slopes. This procedure
is known as the modified Euler’s method.
The modified Euler’s method consists of the following steps.
(i) Predictor step: By using the derivative at the beginning of the step,
the value at the end of the step (t1 = t0 + t) is predicted as
p dx
x1 = x0 + ◊ Dt (5.48)
dt x = x0
p
(ii) Corrector step: By using the predicted value of x1 , the derivative
at the end of the step is computed and the value at the end of the
step (t1 = t0 + t) is predicted as
dx p
= f ( x1 , t1 ) (5.49)
dt p
x = x1
Power System Stability 283
Then, the average value of the two derivatives is used to find the
corrected value.

Ê dx dx ˆ
Á dt + ˜
x = x0 dt x = x1
p
x1 = x0 + Á ˜ ◊ Dt
c
(5.50)
Ë 2 ¯

Similarly,

Ê dx dx ˆ
Á dt + ˜
x = x1 dt p
x = x2
x2 = x1 + Á ˜ ◊ Dt
c
(5.51)
Ë 2 ¯

Ê dx dx ˆ
Á dt + ˜
x = xi dt p
x = xi +1
= xi + Á ˜ ◊ Dt
c
xi + (5.52)
Ë 2 ¯
This process can be repeated until the successive steps converge with
the desired accuracy.

5.8.5 Runge–Kutta (R–K) Method


The Runge–Kutta method approximates the Taylor series solution; however,
unlike the formal Taylor series solution, the R–K method does not require
explicit evaluation of derivatives higher than the first. The effects of the higher
derivatives are included by several evaluations of the first derivative depending
on the number of terms effectively retained in the Taylor series, we have R–K
method of different orders.
Let us consider the first order differential equation
dx
= f ( x, t ) (5.53)
dt
Assume the initial condition is x0, t0.
Second order R–K method
The second order R–K formula for the value of x at t = t0 + t is given by
k1 + k2
x1 = x0 + Dx = x0 + (5.54)
2
where
k1 = (slope at the beginning of time step) t = f (x0, t0)t
k2 = (first approximation to slope at midstep) t = f (x0 + k1, t0 + t)t
284 Electrical Power Systems: Analysis, Security and Deregulation

Fourth order R-K method


The general formula giving the value of x for the (n + 1)th step
1
xi +1 = xi + Dx = xi + (k + 2k2 + 2k3 + k4 ) (5.55)
6 1
where
k1 = (slope at the beginning of time step) t = f (xi, ti)t
Ê k Dt ˆ
k2 = (first approximation to slope at midstep) t = f Á xi + 1 , ti + ˜ Dt
Ë 2 2¯
Ê k Dt ˆ
k3 = (second approximation to slope at midstep) t = f Á xi + 2 , ti + ˜ Dt
Ë 2 2¯
k4 = (slope at the end of step) t = f (xi + k3, ti + t)t

Thus x is the incremental value of x given by the weighted average of estimated


based on slopes at the beginning of the mid-point and end of the time step.
EXAMPLE 5.18 A three-phase fault occurs at the point F as shown in the
figure is cleared by isolating the faulted circuit simultaneously from both
ends. Generator is delivering 0.8 p.u. power at 0.8 p.f lagging. The fault is
cleared in 0.1 s. Obtain the numeric solution of the swing equation up to
0.15 s using the (a) modified Euler’s method and (b) Runge–Kutta method
with step size of t = 0.05 s. Take H = 5 MJ/MVA.

Solution:

(i) To draw the reactance diagram


Et = 1.00°
P = 0.8 p.u.; cosq = 0.8; q = cos–1(0.8) = 36.86°= 0.644 rad
Q = P tanq = 0.6
P - jQ 0.8 - j 0.6
I= = = 0.8 - j 0.6
Et 1.0
Power System Stability 285
E = Et + jXd I = 1.0 + j0.2(0.8 – j0.6) = 1.130.142°
V = Et – j(Xtrans + Xt.line)I = 1.13 + j(0.1 + 0.2) (0.8 – j0.6)
= 0.854–0.285°
(ii) Prefault condition From the reactance diagram,
Ê 0.4 ¥ 0.4 ˆ
X T = X I = 0.2 + 0.1 + Á = 0.5
Ë 0.4 + 0.4 ˜¯
| E ¢| |V | 1.13 ¥ 0.854
Pe1 = sin d = sin d = 1.93 sin d
Xt 0.5

(iii) During the fault condition (fault at the middle of the transmission
line 2)
Using delta to star conversion, the circuit becomes

From the reactance diagram, XT = XII = 1.3


| E ¢| |V | 1.13 ¥ 0.854
Pe 2 = sin d = sin d = 0.742 sin d
X II 1.3
286 Electrical Power Systems: Analysis, Security and Deregulation

(iv) Postfault condition

From the reactance diagram, XT = XIII = j0.7


| E ¢| |V | 1.13 ¥ 0.854
Pe3 = sin d = sin d = 1.378 sin d
X III 0.7

(a) Modified Euler’s method


During fault,
Pe = 0.742 sin d
Pe0 = Pm 0 = 0.8 = 1.93 sin d 0
Ê 0.8 ˆ
d 0 = sin -1 Á = 24.48∞ = 0.427 rad
Ë 1.93 ˜¯
w 0 = 2p f = 2p ¥ 50 = 314.159
Dt = 0.05 s
Iteration 1: Beginning of the first step at, t = 0
dd
= Dw 0 = w 0 - 2p f = 314.159 - 314.159 = 0
dt Dw 0

d Dw pf p ¥ 50
= ( Pm - Pe (d 0 )) = (0.8 - 0.742 sin 0.427) = 15.479
dt d0 H 5
End of the first step, t = 0.05
Predicted values are
p dd
d 0.05 = d 0 + ¥ Dt = 0.427 + 0 ¥ 0.05 = 0.427 rad
dt Dw 0

p d Dw
Dw 0.05 = Dw 0 + ¥ Dt = 0 + 15.479 ¥ 0.05 = 0.774 rad/s
dt d0

Derivatives at the end of t = 0.05


dd p
t = 0.05, = Dw 0.05 = 0.774 rad/s
dt p
Dw 0.05

d Dw pf p p ¥ 50
= ( Pm - Pe (d 0.05 )) = (0.8 - 0.742 sin 0.427) = 15.479 rad/s
dt p
d 0.05 H 5
Power System Stability 287
The corrected values are
È dd ˘ dd
Í dt +
˙ dt
˙ ¥ Dt = 0.427 + È 0 + 0.774 ˘ ¥ 0.05 = 0.446 rad
p
d 0.05 = d 0 + Í
c Dw 0 Dw 0.05
ÍÎ 2 ˙˚ ÍÎ 2 ˙˚
È d Dw d Dw ˘
Í dt +
c Í d0 dt d 0.05 ˙˙
p
È15.479 + 15.479 ˘
Dw 0.05 = Dw 0 + ¥ Dt = 0 + Í ˙˚ ¥ 0.05
ÍÎ 2 ˙
˚ Î 2
= 0.774 rad/s
Iteration 2: Beginning of the first step at t = 0.05,
dd c
= Dw 0.05 = 0.774
dt c
Dw 0.05

d Dw pf c p ¥ 50
= ( Pm - Pe (d 0.05 )) = (0.8 - 0.742 sin 0.446) = 15.077
dt c
d 0.05 H 5
End of the second step, t = 0.1
Predicted values are
p c dd
d 0.1 = d 0.05 + ¥ Dt = 0.446 + 0.774 ¥ 0.05 = 0.485 rad
dt c
Dw 0.05

p c d Dw
Dw 0.1 = Dw 0.05 + ¥ Dt = 0.774 + 15.077 ¥ 0.05 = 1.5278 rad/s
dt c
d 0.05

Derivatives at the end of t = 0.1,
dd p
= Dw 0.1 = 1.5278 rad/s
dt p
Dw 0.1

d Dw pf p p ¥ 50
= ( Pm - Pe (d 0.1 )) = (0.8 - 0.742 sin 0.485) = 14.265 rad/s
dt p
d 0.1 H 5
The corrected values are
È dd dd ˘
Í dt +
dt p ˙
c
È 0.774 + 1.5278 ˘
= d 0.05 + Í
c c Dw 0.05 Dw 0.1 ˙
d 0.1 ¥ Dt = 0.446 + Í ˙˚ ¥ 0.05
ÍÎ 2 ˙˚ Î 2
= 0.5035 rad
È d Dw d Dw ˘
Í dt c + dt p ˙
= Dw 0.05 + Í
c c d 0.05 d 0.1 ˙
Dw 0.1 ¥ Dt
ÍÎ 2 ˙˚
È15.077 + 14.265 ˘
= 0.774 + Í
Î 2 ˙˚ ¥ 0.05 = 1.5076 rad/s
288 Electrical Power Systems: Analysis, Security and Deregulation

The fault is cleared at t = 0.1 s and the postfault condition is given by


Pe = 1.378 sin d
Iteration 3: Beginning of the first step at 0.1 s
dd c
= Dw 0.1 = 1.5076
dt c
Dw 0.1

d Dw pf c p ¥ 50
= ( Pm - Pe (d 0.1 )) = (0.8 - 1.378 sin 0.5035) = 4.245
dt c
d 0.1 H 5

End of the third step, t = 0.15 s


Predicted values are
p c dd
d 0.15 = d 0.1 + ¥ Dt = 0.5035 + 1.5076 ¥ 0.05 = 0.579 rad
dt c
Dw 0.1

p c d Dw
Dw 0.15 = Dw 0.1 + ¥ Dt = 1.5076 + 4.245 ¥ 0.05 = 1.7197 rad/s
dt c
d 0.1

Derivatives at the end of t = 0.15 s,


dd p
= Dw 0.15 = 1.7197 rad/s
dt p
Dw 0.15

d Dw pf p p ¥ 50
= ( Pm - Pe (d 0.15 )) = (0.8 - 1.378 sin 0.579) = 1.444 rad/s
dt p
d 0.15 H 5
The corrected values are
È dd dd ˘
Í dt c + dt p ˙
c
d 0.15
c
= d 0.1 + Í Dw 0.1 Dw 0.15 ˙
¥ Dt
ÍÎ 2 ˙˚
È1.5076 + 1.7197 ˘
= 0.5035 + Í ˙˚ ¥ 0.05 = 0.584 rad
Î 2
È d Dw d Dw ˘
Í dt c + dt p ˙
= Dw 0.1 + Í
c c d 0.1 d 0.15 ˙
Dw 0.15 ¥ Dt
ÍÎ 2 ˙˚
È 4.245 + 1.444 ˘
= 1.5076 + Í
Î 2 ˙˚ ¥ 0.05 = 1.6498 rad/s

(b) Runge–Kutta method
Pe = 0.742 sin d
Pe0 = Pm0 = 0.8 = 1.93 sind0
Power System Stability 289

-1 Ê 0.8 ˆ
d 0 = sin Á = 24.48∞ = 0.427 rad
Ë 1.93 ˜¯
w 0 = 2p f = 2p ¥ 50 = 314.159


Dtt = 0.05 s
Fourth order method
Iteration 1: at t = 0,
Ist estimate: k1 = w0  t = 0  0.05 = 0
pf p ¥ 50
l1 = H ( Pm - Pe (d 0 )) Dt = 5 (0.8 - 0.742 sin 0.427) ¥ 0.05 = 0.774

Ê l ˆ Ê 0.774 ˆ
IInd estimate: k2 = Á Dw 0 + 1 ˜ ¥ Dt = Á 0 + ˜ ¥ 0.05 = 0.0194
Ë 2¯ Ë 2 ¯
pf Ê Ê k1 ˆ ˆ p ¥ 50 Ê Ê 0ˆˆ
l2 = ÁË Pm - Pe ÁË d 0 + ˜¯ ˜¯ Dt = ÁË 0.8 - 0.742 sin ÁË 0.427 + ˜ ˜ ¥ 0.05
H 2 5 2¯¯
= 0.774

IIIrd estimate: Ê l ˆ Ê 0.774 ˆ


k3 = Á Dw 0 + 2 ˜ ¥ Dt = Á 0 + ˜ ¥ 0.05 = 0.0194
Ë 2¯ Ë 2 ¯
pf Ê Ê k ˆˆ p ¥ 50 Ê Ê 0.0194 ˆ ˆ
l3 =Á Pm - Pe Ád 0 + 2 ˜ ˜ Dt = Á 0.8 - 0.742 sin Á 0.427 + ˜ ˜ ¥ 0.05
H Ë Ë 2 ¯ ¯ 5 Ë Ë 2 ¯¯
= 0.764
IVth estimate: k4 = (w0 + l3) × t = (0 + 0.764) × 0.05 = 0.038
pf p ¥ 50
l4 = ( Pm - Pe (d 0 + k3 )) Dt = (0.8 - 0.742 sin (0.427 + 0.0194)) ¥ 0.05
H 5
= 0.7535
Final estimate:
1
d 0.05 = d 0 + (k1 + 2k2 + 2k3 + k4 )
6
1
= 0.427 + (0 + 2 ¥ 0.0194 + 2 ¥ 0.0194 + 0.038)
6
= 0.446 rad
1
Dw 0.05 = Dw 0 + (l1 + 2l2 + 2l3 + l4 )
6
1
= 0 + (0.774 + 2 ¥ 0.774 + 2 ¥ 0.764 + 0.7535)
6
= 0.767 rad/s

290 Electrical Power Systems: Analysis, Security and Deregulation

Iteration 2: at t = 0.05s,
Ist estimate: k1 = w0.05  t = 0.767  0.05 = 0.0384
pf p ¥ 50
l = ( Pm - Pe (d 0.05 )) Dt = (0.8 - 0.742 sin 0.446) ¥ 0.05 = 0.754
1 H 5
Ê l ˆ Ê 0.754 ˆ
IInd estimate: k2 = Á Dw 0.05 + 1 ˜ ¥ Dt = Á 0.767 + ˜ ¥ 0.05 = 0.0572
Ë 2¯ Ë 2 ¯
pf Ê Ê k ˆˆ
l2 = Á Pm - Pe Ád 0.05 + 1 ˜ ˜ Dt
H Ë Ë 2 ¯¯
p ¥ 50 Ê Ê 0.0384 ˆ ˆ
= Á 0.8 - 0.742 sin Á 0.446 + ˜ ˜ ¥ 0.05
5 Ë Ë 2 ¯¯
= 0.734

IIIrd estimate: Ê l ˆ Ê 0.734 ˆ


k3 = Á Dw 0.05 + 2 ˜ ¥ Dt = Á 0.767 + ˜ ¥ 0.05 = 0.0567
Ë 2 ¯ Ë 2 ¯
pf Ê Ê k ˆˆ
l3 = Á Pm - Pe Ád 0.05 + 2 ˜ ˜ Dt
H Ë Ë 2 ¯¯
p ¥ 50 Ê Ê 0.0 572ˆ ˆ
= Á 0.8 - 0.742 sin Á 0.446 + ˜ ˜ ¥ 0.05
5 Ë Ë 2 ¯¯
= 0.724
IVth estimate: k4 = (w0.05 + l3)  t = (0.767 + 0.724)  0.05 = 0.0746
pf
l4 = (Pm - Pe (d 0.05 + k3 )) Dt
H
p ¥ 50
= (0.8 - 0.742 sin (0.446 + 0.0567)) ¥ 0.05
5
= 0. 695
Final estimate:
1
d 0.1 = d 0.05 + (k1 + 2k2 + 2k3 + k4 )
6
1
= 0.446 + (0.0384 + 2 ¥ 0.0572 + 2 ¥ 0.0567 + 0.0746)
6
= 0.503 rad
1
Dw 0.1 = Dw 0.05 + (l1 + 2l2 + 2l3 + l4 )
6
1
= 0.767 + (0.754 + 2 ¥ 0.734 + 2 ¥ 0.724 + 0.695)
6
= 1.495 rad/s
Fault is cleared at t = 0.1
Postfault condition, Pe = 1.378 sin d
Power System Stability 291
Iteration 3:
Ist estimate: k1 = w0.1  t = 1.495  0.05 = 0.0748
pf p ¥ 50
l1 = ( Pm - Pe (d 0.1 )) Dt = (0.8 - 1.378 sin 0.503) ¥ 0.05 = 0.213
H 5
Ê l ˆ Ê 0.213 ˆ
IInd estimate: k2 = Á Dw 0.1 + 1 ˜ ¥ Dt = Á1.495 + ˜ ¥ 0.05 = 0.08
Ë 2 ¯ Ë 2 ¯
pf Ê Ê k1 ˆ ˆ
l2 = ÁË Pm - Pe ÁË d 0.1 + ˜¯ ˜¯ Dt
H 2
p ¥ 50 Ê Ê 0.748ˆ ˆ
= ÁË 0.8 - 1.378 sin ÁË 0.503 + ˜ ˜ ¥ 0.05
5 2 ¯¯
= 0.143
Ê l ˆ Ê 0.143 ˆ
IIIrd estimate: k3 = Á Dw 0.1 + 2 ˜ ¥ Dt = Á1.495 + ˜ ¥ 0.05 = 0.0783
Ë 2 ¯ Ë 2 ¯
pf Ê Ê k ˆˆ
l3 = Á Pm - Pe Ád 0.1 + 2 ˜˜ Dt
H Ë Ë 2 ¯¯
p ¥ 50 Ê Ê 0.0 8ˆ ˆ
= Á 0.8 - 1.378sin Á 0. 503 + ˜ ˜ ¥ 0.05
5 Ë Ë 2 ¯¯
= 0. 1382
IVth estimate: k4 = (w0.1 + l3)  t = (1.495 + 0.1382)  0.05 = 0.0817
pf
l4 = (Pm - Pe (d 0.1 + k3 )) Dt
H
p ¥ 50
= (0.8 - 1.378 sin (0. 503 + 0.0783)) ¥ 0.05
5
= 0. 068
Final estimate:
1
d 0.15 = d 0.1 + (k1 + 2k2 + 2k3 + k4 )
6
1
= 0.503 + (0.0748 + 2 ¥ 0.08 + 2 ¥ 0.0783 + 0.0817)
6
= 0.582 rad
1
Dw 0.15 = Dw 0.1 + (l1 + 2l2 + 2l3 + l4 )
6
1
= 1.495 + (0.213 + 2 ¥ 0.143 + 2 ¥ 0.1382 + 0.068)
6
= 1.636 rad/s
292 Electrical Power Systems: Analysis, Security and Deregulation

5.9 Multimachine Transient Stability


Transient stability analysis has recently become a major issue in the operation
of power system due to the increasing stress on power system networks. This
problem requires evaluation of a power system’s ability to withstand disturbances
while maintaining the quality of service. Many different techniques have been
proposed for transient stability analysis in power system, especially for a
multi-machine system.
Multi-machine equations can be written similar to the one-machine system
connected to the infinite bus. In order to reduce the complexity of the transient
stability analysis, similar simplifying assumptions are made as follows.
• Each synchronous machine is represented by a constant voltage source
behind the direct axis transient reactance. This representation neglects
the effect of saliency and assumes constant flux linkages.
• The actions of the governor are neglected and the input powers are
assumed to remain constant during the entire period of simulation.
• Using the prefault bus voltages, all loads are converted to equivalent
admittances to ground and are assumed to remain constant.
• Damping or asynchronous powers are ignored.
• The mechanical rotor angle of each machine coincides with the angle
of the voltage behind the machine reactance.
• Machines belonging to the same station swing together and are said
to be coherent. A group of coherent machines is represented by one
equivalent machine.

5.9.1 Mathematical Model of Multimachine Transient


Stability Analysis
The first step in the transient stability analysis is to solve the initial load
flow and to determine the initial bus voltage magnitudes and phase angles.
The machine currents prior to disturbance are calculated from,

Si* Pi - jQi
Ii = = i = 1, 2, ..., m
Vi* Vi* (5.56)
where m is the number of generators
Vi is the terminal voltage of the ith generator
Pi and Qi are the generators of real and reactive powers.
All unknown values are determined from the initial power flow solution.
The generator armature resistances are usually neglected and the voltages
behind the transient reactance are then obtained as
Ei¢ = Vi + jX d¢ I i (5.57)
Power System Stability 293
Next, all loads are converted to equivalent admittances by using the relation

Si* Pi - jQi
yi 0 = 2
= (5.58)
|Vi | |Vi |2
To include voltages behind the transient reactance, m buses are added to the n
bus power system network. The equivalent network, with all loads converted
to admittances is shown in Figure 5.30.
Nodes n + 1, n + 2, ..., n + m are the internal machine buses, i.e. the
buses behind the transient reactances. The node voltage equation, with node
0 as reference for this network, is

Figure 5.30 Power system representation for transient stability analysis.

È I1 ˘ È Y11  Yn1 Y1( n +1)  Y1( n + m ) ˘ È V1 ˘


Í  ˙ Í       ˙ Í  ˙
Í ˙ Í ˙ Í ˙
Í I n ˙ Í Y1n  Ynm Yn ( n +1)  Yn ( n + m ) ˙ Í Vn ˙
Í ˙=Í ˙ Í ˙
Í I n +1 ˙ Í Y( n +1)1  Y( n +1) n Y( n +1) ( n +1)  Y( n +1) ( n + m ) ˙ Í En¢ +1 ˙
Í  ˙ Í       ˙ Í  ˙
Í ˙ Í ˙ Í ˙
ÍÎ I n + m ˙˚ ÍÎY( n + m )1  Y( n + m ) n Y( n + m ) ( n +1)  Y( n + m ) ( n + m ) ˙˚ ÍÎ En¢ + m ˙˚
(5.59)
or
Ibus = YbusVbus (5.60)
where Ibus is the vector of the injected bus currents
Vbus is the vector of bus voltages measured from the reference node.
The diagonal elements of the bus admittance matrix are the sum of
admittances connected to it, and the off-diagonal elements are equal to the
negative of the admittance between the nodes. The reference is that additional
nodes are added to include the machine voltages behind transient reactances.
Also, the diagonal elements are modified to include the load admittances.
294 Electrical Power Systems: Analysis, Security and Deregulation

To simplify the analysis, all nodes other than the generator internal nodes
are eliminated using Kron’s reduction formula. To eliminate the load buses,
the bus admittance matrix in Eq. (5.59) is partitioned such that the n buses to
be removed are represented in the upper n rows. Since no current enters or
leaves the load buses, currents in the n rows is zero. The generator current is
denoted by the vector Im and the generator and load voltages are represented
by the vectors Em and Vn, respectively. Then, Eq. (5.59), in terms of sub
matrices becomes

È 0 ˘ È Ynn Ynm ˘ È Vn ˘
ÍI ˙ = Í t ˙ ÍE¢ ˙ (5.61)
Î m ˚ ÎYnm Ymm ˚ Î m˚

The voltage vector Vn may be eliminated upon substitution as follows.


0 = YnnVn + Ynm Em¢ (5.62)
t
Im = YnmVn + Ymm Em¢ (5.63)
From Eq. (5.62)
-1
Vn = - Ynn + Ynm Em¢ (5.64)
Now substituting into Eq. (5.63)
t -1 red
I m = [Ymm - YnmYnn Ynm ] Em¢ = Ybus Em¢ (5.65)
The reduced admittance matrix is
red t -1
Ybus = Ymm - YnmYnn Ynm (5.66)
The reduced bus admittance matrix has the dimensions (m  m), where m is
the number of generators. The electrical power output of each machine can
now be expressed in terms of the machine’s internal voltages
Sei* = E'i*Ii
or
Pei = Re(E'i*Ii) (5.67)
where
m
Ii = Â E ¢Y
j =1
j ij (5.68)

Expressing voltages and admittances in the polar form,

Ei¢ = |Ei¢|–d i and Yij = |Yij |–qij


Substituting Ii in Eq. (5.67), we get
m
Pei = Â |E ¢| |E ¢ | |Y | cos(q
j =1
i j ij ij - di + d j ) (5.69)
Power System Stability 295
The above equation is the same as the power flow equation. Prior to disturbance,
there is equilibrium between the mechanical power input and the electrical
power output, and we have
m
Pmi = Â |E ¢| |E ¢ | |Y | cos(q
j =1
i j ij ij - di + d j ) (5.70)

The classical transient stability study is based on the application of a three-


phase fault. A solid three-phase fault at bus k in the network results in Vk = 0.
This is simulated by removing the kth row and column from the prefault bus
admittance matrix. The new bus admittance matrix is reduced by eliminating
all nodes except the internal generator nodes. The generator excitation voltages
during the fault and the postfault modes are assumed to remain constant.
The electrical power of the ith generator in terms of the new reduced bus
admittance matrices are obtained from Eq. (5.69). The swing equation with
damping neglected, for machine i becomes
2 m
Hi d di

p f 0 dt 2
= Pmi - Â |E ¢| | E ¢ | |Y | cos(q
j =1
i j ij ij - di + d j ) (5.71)

where, Yij are the elements of the faulted reduced bus admittance matrix,
Hi is the inertia constant of machine i expressed on the common MVA base SB.
If HGi is the inertia constant of machine i expressed on the machine rated
MVA SGi, then Hi is given by
SGi
Hi = H Gi (5.72)
Scb
Showing the electrical power of the ith generator by Pef and transforming
Eq. (5.71) into state variable mode yield
dd i
= wi (5.73)
dt
d Dw i p f0 f
= ( Pmi - Pei ) (5.74)
dt Hi
In the transient stability analysis problem, we have two state equations for
each generator. When the fault is cleared, which may involve the removal of
the faulty line, the bus admittance matrix is recomputed to reflect the change
in the network. Next the postfault reduced bus admittance matrix is evaluated
pf
and the postfault electrical power of the ith generator shown by Pi readily
pf
determined. Using the postfault power Pi , the simulation is continued to
determine the system stability, until the plots reveal a definite trend as to
stability or instability. Usually the slack generator is selected as the reference
machines are plotted. Usually, the solution is carried out for two swings to show
296 Electrical Power Systems: Analysis, Security and Deregulation

that the second swing is not greater than the first one. If the angle differences
do not increase, the system is stable. If any of the angle differences increase
indefinitely, the system is unstable. The flow chart of transient stability analysis
for a multimachine power is given in Figure 5.31.

Figure 5.31 Flow chart of transient stability analysis for a multimachine power.

5.10 Factors Influencing Transient Stability


The factors that affect the transient stability are given below.
(i) The generator inertia. The higher the inertia, the lower is the rate
of change in angle. This reduces the kinetic energy gained during
fault.
(ii) The generator reactance. A lower reactance increases peak power and
reduces initial rotor angle.
(iii) The generator internal voltage magnitude (E ). This depends on the
field excitation.
(iv) How heavily the generator is loaded.
(v) The generator output during the fault. This depends upon the fault
location and type.
(vi) The fault clearing time.
(vii) The postfault transmission system reactance.
Power System Stability 297

5.11 Techniques for Transient Stability


Improvement
(i) Reduction in system transfer reactance
Reducing the series reactance by using the series capacitor is normally
economical for lines of length more than 320 km. For lines of length less
than 320 km, the objective is achieved by running parallel lines. When parallel
lines are used, instead of a single line, some power can be transferred over
the healthy line even during a three-phase fault on one of the lines, unless of
course when a fault takes place at the paralleling bus when no power can be
transferred out of the parallel lines. For other types of faults on more than one
line, more power can be transferred during the fault, if there are two lines in
parallel, power can be transferred over a single faulted line.
The effect of reducing the series reactance is to increase Pm, which
therefore, increases the transient stability of a system.
(ii) Increase of system voltage
From the equation Pmax = EGEM /XT, it is clear that an increase in system voltage
results in higher values of power Pmax that can be transferred between the nodes.
Since shaft power Pe = Pmax sin d0, with higher values of Pmax, d0 is reduced
and, therefore, the difference between the critical clearing angle and the initial
angle d0 is increased. Therefore, increasing Pmax allows the machine to rotate
through large angle before it reaches the critical clearing angle, which results
in greater critical clearing time and the probability of maintaining stability.
(iii) Use of high speed reclosing breakers
The quicker a breaker operates, the faster the fault is removed from the system
and the better is the tendency of the system to restore to normal operating
conditions. The use of high-speed breakers has materially improved the transient
stability of the power systems and does not require any other methods for
the purpose. The use of reclosing type circuit breakers plays a vital role in
improving the transient stability limit.
(iv) HVDC links
A dc link is asynchronous, i.e. the two ac systems at either end do not have
to be controlled in phase or even be at exactly the same frequency as they
do for an ac link. There is no risk of a fault in one system causing loss of
stability in the other system.
(v) Braking resistors
For improving stability, when large load is suddenly lost, a resistive load
called a braking resistor is connected at or near the generator bus. This load
compensates for at least some of the reduction of load on the generators and
so reduces the acceleration. During a fault, the resistors are applied to the
terminals of the generators through circuit breakers by means of elaborate
298 Electrical Power Systems: Analysis, Security and Deregulation

control schemes. The control scheme determines the amount of resistance to


be applied and its duration.
(vi) Short circuit current limiters
These may be used in long transmission lines to modify favourably the transfer
impedance during the fault conditions so that the voltage profile of the system
is somewhat improved, thereby raising the system load level during the fault.
(vii) Turbine fast valving or bypass valving
Another recent method of improving the stability of a unit is to decrease the
mechanical input power to the turbine. This can be accomplished by means
of fast valving, where the difference between the mechanical input and the
reduced electrical output of a generator under a fault, as sensed by a control
scheme, initiates the closing of a turbine valve to reduce the power input.
(viii) Full load rejection technique
In places where stability is difficult to maintain, the normal procedure is to
automatically trip the unit off the line. This, however, causes several hours
of delay before the unit can be put back into operation. The loss of a major
unit for this length of time can seriously jeopardize the remaining system.
To remedy these situations, a full load rejection scheme could be utilized after
the unit is separated from the system. To do this, the unit has to be equipped
with a large steam bypass system. After the system has recovered from the
shock caused by the fault, the unit could be resynchronized and reloaded.
The main disadvantage of this method is the extra cost of a large bypass system.

Review Questions
Part-A
1. What is the power system stability?
2. How is the power system stability classified?
3. What is the rotor angle stability?
4. What is the steady state stability?
5. What is the steady state stability limit?
6. What is the transient stability?
7. What is the transient stability limit?
8. What is the dynamic stability?
9. What is the voltage stability?
10. State the causes of voltage instability.
11. Write the power angle equation and draw the power angle curve.
12. Write the expression for the maximum power transfer.
13. Write the swing equation for a SMIB (single machine connected to an
infinite bus bar) system.
Power System Stability 299
14. Define the swing curve.
15. In a three machine system having ratings G1, G2 and G3 and inertia
constants M1, M2 and M3 what are the inertia constants M and H of
the equivalent system.
16. State the assumptions made in stability studies.
17. State equal area criterion.
18. Define the critical clearing angle.
19. List the methods of improving the transient stability limit of a power
system.
20. What are the numerical integration methods of power system stability?

Part-B
1. A 400 MVA synchronous machine has H1 = 4.6 MJ/MVA and a
1200 MVA machine H2 = 3.0 MJ/MVA. Two machines operate in
parallel in a power plant. Find out Heq relative to a 100 MVA base.
2. A 100 MVA, two pole, 50 Hz generator has moment of inertia
40  103 kg · m2. What is the energy stored in the rotor at the rated
speed? What is the corresponding angular momentum? Determine the
inertia constant H.
3. The sending end and the receiving end voltages of a three-phase
transmission line at a 200 MW load are equal at 230 kV. The per phase
line impedance is j14 . Calculate the maximum steady state power
that can be transmitted over the line.
4. A single line diagram of a system is shown in the figure below.
All the values are in per unit on a common base. The power delivered
into bus 2 is 1.0 p.u. at 0.80 power factor lagging. Obtain the power
angle equation and the swing equation for the system. Neglect all losses.

  

5. A 50 Hz synchronous generator capable of supplying 400 MW of power


is connected to a larger power system and is delivering 80 MW when
a three-phase fault occurs at its terminals, determine (a) the time in
which the fault must be cleared if the maximum power angle is to be
–85°, assume H = 7 MJ/MVA on a 100 MVA base and (b) the critical
clearing angle.
6. A synchronous generator is connected to a large power system and
supplying 0.45 p.u. MW of its maximum power capacity. A three-
300 Electrical Power Systems: Analysis, Security and Deregulation

phase fault occurs and the effective terminal voltage of the generator
becomes 25% of its value before the fault. When the fault is cleared, the
generator is delivering 70% of the original maximum value. Determine
the critical clearing angle.
7. Determine the critical clearing angle of the power system as shown in
the figure below for a three-phase fault at the point F. The generator
is supplying 1.0 p.u. MW power under the prefault condition.

You might also like