Chapter 5
Chapter 5
Chapter 5
Introduction Stability of a power system is its ability to return to normal or stable operating conditions after having been subjected to some form of disturbance. Conversely, instability means a condition denoting loss of synchronism or falling out of step. Though stability of a power system is a single phenomenon, for the purpose of analysis, it is classified as Steady State Analysis and Transient Stability. Increase in load is a kind of disturbance. If increase in loading takes place gradually and in small steps and the system withstands this change and performs satisfactorily, then the system is said to be in STADY STATE STABILITY. Thus the study of steady steady state stability is basically concerned with the determination of upper limit of machines loading before losing synchronism, provided the loading is increased gradually at a slow rate.
In practice, load change may not be gradual. Further, there may be sudden disturbances due to i) ii) iii) iv) Sudden change of load Switching operation Loss of generation Fault
Following such sudden disturbances in the power system, rotor speeds, rotor angular differences and power transfer undergo fast changes whose magnitudes are dependent upon the severity of disturbances, For a large disturbance, changes in angular differences may be so large as to cause the machine to fall out of step. This type of instability is known as TRANSIENT INSTABILITY. Transient stability is a fast phenomenon, usually occurring within one second for a generator close to the cause of disturbance.
Short circuit is a severe type of disturbance. During a fault, electrical powers from the nearby generators are reduced drastically, while powers from remote generators are scarily affected. In some cases, the system may be stable even with sustained fault; whereas in other cases system will be stable only if the fault is cleared with sufficient rapidity. Whether the system is stable on the occurrence of a fault depends not only on the system itself, but also on the type of fault, location of fault, clearing time and the method of clearing. Transient stability limit is almost always lower than the steady state limit and hence it is much important. Transient stability limit depends on the type of disturbance, location and magnitude of disturbance. Review of mechanics Transient stability analysis involves some mechanical properties of the machines in the system. After every disturbance, the machines must adjust the relative angles of their rotors to meet the condition of the power transfer involved. The problem is mechanical as well as electrical.
(1)
where I is the Moment of Inertia in Mega Joules sec.2 / elec. deg.2 is the angular velocity in elec. deg. / sec. Angular Momentum M = I ; Then from eqn. (1), K.E. can be written as K.E. =
1 M Mega Joules 2
(2)
The angular momentum M depends on the size of the machine as well as on its type. The Inertia constant H is defined as the Mega Joules of stored energy of the machine at synchronous speed per MVA of the machine. When so defined, the relation between the Angular Momentum M and the Inertia constant H can be derived as follows.
Let G be the rating of the machine in MVA. Then Stored energy = G H MJ Further, K.E. =
1 M MJ 2
(3) =
1 M (2 f) MJ = M x f MJ 2
(4)
(5)
If the power is expressed in per unit, then G = 1.0 per unit and hence M=
H f
(6)
While the angular momentum M depend on the size of the machine as well as on its type, inertia constant H does not vary very much with the size of the machine, The quantity H has a relatively narrow range of values for each class of machine.
Swing equation The differential equation that relates the angular momentum M, the acceleration power Pa and the rotor angle is known as SWING EQUATION. Solution of swing equation will show how the rotor angle changes with respect to time following a disturbance. The plot of Vs t is called the SWING CURVE. Once the swing curve is known, the stability of the system can be assessed. The flow of mechanical and electrical power in a generator and motor are shown in Fig. 1.
Te
Generator
Pe Ts
Motor
Pe
Pm Ts
(a)
Pm Te
Fig. 1
(b)
Consider the generator shown in Fig. 1(a). It receives mechanical power P m at the shaft torque T s and the angular speed via. shaft from the prime -mover. It delivers electrical power Pe to the power system network via. the bus bars. The generator develops electromechanical torque T e in opposition to the shaft torque Ts. At steady state, T s = T e.
Assuming that the windage and the friction torque are negligible, the accelerating torque acting on the rotor is given by Ta = Ts Te Multiplying by on both sides, we get Pa = Ps Pe In case of motor Ta = Te Ts Pa = Pe Ps In general, the accelerating power is given by Pa = Input Power Output Power (11) (9) (10) (8) (7)
d2 d = = angular acceleration = dt dt 2
We know that
d2 Pa = T a = I = M = M dt 2
(12)
Consider an object moving at a linear speed of v s v. It is required to find its displacement at any time t. For this purpose, introduce another object moving with a constant speed of vs. Then at any time t, the diaplacement of the first object is given by x = vs t + d where d is the displacement of the first object with respect to the second object as shown in Fig. 2.
x vs + v vs d
Fig. 2
Similarly in the case of angular movement, the angular displacement , at any time t is given by
= s t +
(13)
where is the angular displacement of the rotor from synchronously rotating reference axis, rotating at synchronous speed s. The angle is also called as LOAD ANGLE or TORQUE ANGLE. In view of eqn.(13)
d d = s + dt dt
(14)
d2 d2 = dt 2 dt 2
(15)
(16)
(17)
Swing curve, which is the plot of torque angle vs time t, can be obtained by solving the swing equation. Two typical swing curves are shown in Fig. 3.
.
s
d 0 dt
0
t (a)
0
t (b)
Fig. 3 Swing curves are used to determine the stability of the system. If the rotor angle
reaches a maximum and then decreases, then it shows that the system has
transient stability. On the other hand if the rotor angle increases indefinitely, then it shows that the system is unstable.
We are going to study the stability of (1) a generator connected to infinite bus and (2) a synchronous motor drawing power from infinite bus. We know that the complex power is given by P + j Q = V I* i.e. P j Q = V * I Thus real power P = Re {V * I}
E j X I = V i.e. E = V + j X I
E V I
Current I =
E V X
Xd M E I
XT
VjXI = E
V I
V E
Thus E = E
Current I =
1 [ V ( E cos j E sin )] jX
E V X
Notice that the swing equation is second order nonlinear differential equation
Equal area criterion The accelerating power in swing equation will have sine term. Therefore the swing equation is non-linear differential equation and obtaining its solution is not simple. For two machine system and one machine connected to infinite bus bar, it is possible to say whether a system has transient stability or not, without solving the swing equation. Such criteria which decides the stability, makes use of equal area in power angle diagram and hence it is known as EQUAL AREA CRITERION. Thus the principle by which stability under transient conditions is determined without solving the swing equation, but makes use of areas in power angle diagram, is called the EQUAL AREA CRITERION. From the Fig. 3, it is clear that if the rotor angle oscillates, then the system is stable. For to oscillate, it should reach a maximum value and then should decrease. At that point
d = 0. Because of damping inherently present in the dt
system, subsequence oscillations will be smaller and smaller. Thus while changes, if
d = 0, then the stability is ensured. dt
The swing equation for the alternator connected to the infinite bus bars is
d2 M = Ps Pe dt 2
(18)
d , we get dt
i.e.
1 d d d M ( ) 2 (Ps Pe) 2 dt dt dt
(19)
Thus
d d 2 dt 2 (Ps Pe ) ( ) ; dt dt d M
i.e.
d d 2 2 (Ps Pe ) ( ) d dt M
On integration
d ( )2 = dt
2 (Ps Pe ) d M 0
i.e.
d = dt
2 (Ps Pe ) d M 0
(20)
Before the disturbance occurs, 0 was the torque angle. At that time the disturbance occurs,
d is no longer zero and starts changing. dt
d = 0. Once dt
Torque angle will cease to change and the machine will again be operating at synchronous speed after a disturbance, when
d = 0 or when dt
2 (Ps Pe ) d = 0 i.e. M 0
(P
0
Pe ) d = 0
(21)
If there exist a torque angle for which the above is satisfied, then the machine will attain a new operating point and hence it has transient stability. The machine will not remain at rest with respect to infinite bus at the first time when
d = 0. But due to damping present in the system, during subsequent dt
oscillation, maximum value of keeps on decreasing. Therefore, the fact that has momentarily stopped changing may be taken to indicate stability.
Sudden load increase on Synchronous motor Let us consider a synchronous motor connected to an infinite bus bars.
M I V E Pe = E jX V E = V I ( jX )
E V X
Ps
Fig. 4
P c
Ps
d A1
b A2 e Changed output
P0
Fig. 4
The following changes occur when the load is increased suddenly. Point a Initial condition; Input = output = P0; = s; = 0 Due to sudden loading, output = Ps; output > Input; decreases from s; increases from 0. Between a-b Output > Input; Deceleration; decreases; increases. Point b Output = Input; = min which is less than s; = s Since is less than s, continues to increase.
Ps
Between b-c Input > output; Rotating masses start gaining energy; Acceleration; starts increasing from minimum value but still less than s; continues to increase. Point c Input > output; = s; = m; There is acceleration; is going to increase from s; hence is going to decrease from m. Between c-b Input > output; Acceleration; increases and decreases. Point b Input = output; = max ; = s. Since is greater than s, continues to decrease. Between b-a Output > input; Deceleration; starts decreasing from max ; but still greater than s; continues to decrease. Point a = s; = 0; Output > Input; The cycle repeats.
Because of damping present in the system, subsequent oscillations become smaller and smaller and finally b will be the steady state operating point. Interpretation of equal area As discussed earlier (eqn. 21), the condition for stability is
(Ps Pe ) d = 0 i.e.
m
Pe d =
Ps d
From Fig. 4,
P
0
d = area 0 a b c m
Ps d A1
b A2 e Changed output
and
P
0
d = area 0 a d e m
P0 a
Initial output
Fig. 4
Subtracting area 0 a b e m from both sides of above equation, we get A2 = A1. Thus for stability, A2 = A1 (22)
Fig. 5 shows three different cases: The one shown in case a is STABLE. Case b indicates CRITICALLY STABLE while case c falls under UNSTABLE.
P P
A2 A2 A1 A1 A1
P
A2
P0
P0
P0
0 s
0 s
0
Case c
m
Case a
Case b Fig. 5 5
Example 1 A synchronous motor having a steady state stability limit of 200 MW is receiving 50 MW from the infinite bus bars. Find the maximum additional load that can be applied without causing instability.
P
200 PS
A1
A2
50
B A E
- S
Fig. 6 5
0 = sin-1
Further 200 sin S = PS Adding area ABCDEA to both A1 and A2 and equating the resulting areas
200 sin S ( S 0) =
200 sin d
i.e.
( S 0) sin S = cos 0 cos ( S) = cos 0 + cos S i.e. ( S 0.25268) sin S - cos S = 0.9682458 The above equation can be solved by trial and error method. S RHS 0.85 0.8718 0.9 0.9363 0.95 0.9954
Using linear interpolation between second and third points we get S = 0.927 rad. 0.927 rad. = 53.11 deg. Thus PS = 200 sin 53.110 = 159.96 MW Maximum additional load possible = 159.96 50 = 109.96 MW
Opening of one of the parallel lines When a generator is supplying power to an infinite bus over two parallel transmission lines, the opening of one of the lines will result in increase in the equivalent reactance and hence decrease in the maximum power transferred. Because of this, depending upon the initial operating power, the generator may loose synchronism even though the load could be supplied over the remaining line under steady state condition. Consider the system shown in Fig. 7. The power angle diagrams corresponding to stable and unstable conditions are shown in Fig. 8.
G Fig. 7 5
A2
PS
A2 A1
PS
A1
s m
A2 = A1 ;
Stable
A2 < A1 ;
Fig. 8
Unstable
Short circuit occurring in the system Short circuit occurring in the system often causes loss of stability even though the fault may be removed by isolating it from the rest of the system in a relatively short time. A three phase fault at one end of a double circuit line is shown in Fig. 9(a) which can be reduced as shown in Fig. 9(b).
+ + Eg (a)
+ Ei Eg Ig
+ Ei -
(b) Fig. 9 5 generator flows through the fault and It is to be noted that all the current from the
this current Ig lags the generator voltage by 90 0. Thus the real power output of the generator is zero. Normally the input power to the generator remains unaltered. Therefore, if the fault is sustained, the load angle will increase indefinitely because all the input power will be used for acceleration, resulting unstable condition.
When the three phase fault occurring at one end of a double circuit line is disconnected by opening the circuit breakers at both ends of the faulted line, power is again transmitted. If the fault is cleared before the rotor angle reaches a particular value, the system will remain stable; otherwise it will loose stability as shown in Fig. 10.
A2
A2
Input power
PS
A1
A2 A1
A1
0 C m
C < CC
CC
A2 = A1 ;
Stable
A2 = A1 ; Critically stable
Fig. 10
When a three phase fault occurs at some point on a double circuit line, other than on the extreme ends, as shown in Fig. 11(a), there is some finite impedance between the paralleling buses and the fault. Therefore, some power is transmitted during the fault and it may be calculated after reducing the network to a delta connected circuit between the internal voltage of the generator and the infinite bus as shown in Fig. 11(b). Xb
+ + Eg (a) Fig. 11
+ Em Eg -
Xa
Xc
+ Em -
(b)
Eg
Em Xb
sin
(23)
Stable, critically stable and unstable conditions of such systems are shown:
CRITICALLY STABLE
CC
UNSTABLE
A1
Example 2 In the power system shown in Fig. 12, three phase fault occurs at P and the faulty line was opened a little later. Find the power output equations for the pre-fault, during fault and post-fault conditions.
0.16 Xg = 0.28 G P x 0.16 0.24 0.16
E g 1.25 p.u.
'
V 1.0 p.u.
= 0.28
Fig. 12
0.72 + 1.0 -
Power output Pe =
During fault condition: 0.16 0.56 0.16 0.28 + 1.25 0.4 + 1.0 0.28 + 0.16 0.4 0.16 + 0.56
0.28 + 1.25 -
0.08
0.2 + 0.057 -
Power output Pe =
Post-fault condition: 0.16 0.56 1.0 0.28 + 1.25 + 1.0 1.25 1.0
Power output Pe =
Thus power output equations are: Pre-fault During fault Post fault Here Pm1 = 1.736; Pm2 = 0.418; Pm3 = 1.25; Pe = Pm1 sin = 1.736 sin Pe = Pm2 sin = 0.418 sin Pe = Pm3 sin = 1.25 sin
C C
0
m2
CC Fig. 13
Area A1 = PS ( C C 0 )
sin d
= PS C C PS 0 Pm 2 cos C C Pm 2 cos 0
m
(24)
Area A2 =
C C
m3
sin d - PS ( m C C )
= Pm 3 cos C C Pm 3 cos m PS m PS C C
(25)
(24) (25)
Area A2 = Area A1
Pm 3 cos C C Pm 3 cos m PS m PS C C = PS C C PS 0 Pm 2 cos C C Pm 2 cos 0 ( Pm 3 Pm 2 ) cos C C PS ( m 0 ) Pm 2 cos 0 Pm 3 cos m
cos C C
PS ( m 0 ) Pm 2 cos 0 Pm 3 cos m Pm 3 Pm 2
C C cos 1 [
PS ( m 0 ) Pm 2 cos 0 Pm 3 cos m Pm 3 Pm 2
(26)
m sin 1 (
PS ) Pm 3
(27)
Example 3 In the power system described in the previous example, if the generator was delivering 1.0 p.u. just before the fault occurs, calculate C C . Solution Pm1 = 1.736; Pm2 = 0.418; Pm3 = 1.25; PS = 1.0
1.736 sin 0 = 1.0; sin 0 = 0.576; 0 = 0.6139 rad. 1.25 sin ( - m ) = 1.0; sin ( - m ) = 0.8; - m = 0.9273; m = 2.2143 rad.
cos C C
PS ( m 0 ) Pm 2 cos 0 Pm 3 cos m Pm 3 Pm 2
1.0 (2.2143 0.6139) 0.418 cos 0.6139 1.25 cos 2.2143 0.6114 1.25 0.418
STEP BY STEP SOLUTION OF OBTAINING SWING CURVE The equal area criterion of stability is useful in determining whether or not a system will remain stable and in determining the angle through which the machine may be permitted to swing before a fault is cleared. It does not determine directly the length of time permitted before clearing a fault if stability is to be maintained. In order to specify a circuit breaker of proper speed, the engineer must know the CRITICAL CLEARING TIME, which is the time for a machine to swing from its original position to its critical clearing angle. If the Critical Clearing Angle (CCA) is determined by the equal area criterion, then to determine corresponding Critical Clearing Time (CCT), the swing curve for the sustained faulted condition is required. The step by step method of obtaining swing curve, using hand calculation is necessarily simpler than some of the methods recommended for digital computer. In the method for hand calculation, the change in the angular position of the rotor during a short interval of time is computed by making the following assumptions.
1. The accelerating power Pa computed at the beginning of an interval is constant from the middle of the proceeding interval to the middle of the interval considered. 2.
d is constant throughout any interval at the value computed for the dt
middle of the interval. Fig. 14 will help in visualizing the assumptions. The accelerating power is computed for the points enclosed in circles, at the beginning of n-1, n and n+1 th intervals. The step of Pa in the figure results from assumption 1. Similarly (
d d = - s), the excess of angular velocity over the synchronous dt dt
angular velocity is shown as a step curve that is constant throughout the interval, at the value computed at the midpoint. Between the ordinates n -
Pa (n 2) Pa (n 1) Pa (n) t t n-2 (n 1/2) (n 3/2) t t n n-1 n-2 n-2 n-1 Fig. 14 n n n-1 n 3/2 n 1/2 n-1 n t
n th interval n - 1 th interval
Pa ( n 1) M
d2 1 ( Because = Pa ) 2 M dt
(28)
Similarly, change in over any interval = constant angular speed x time duration. Thus (n 1) = (n 3 ) t 2
and
(29)
(n) = (n -
1 ) t 2
Pa ( n 1) 1 3 ) - (n - ) ] t = (t)2 M 2 2
(30)
Therefore
(n) - (n 1) = [ (n Pa ( n 1) M
Thus
(n) = (n 1) +
(t)2
(31)
Thus
(n) = (n 1) +
Pa ( n 1) M
(t)2
(31)
Equation (31) shows that the change in torque angle during a given interval is equal to the change in torque angle during the proceeding interval plus the
( t) 2 accelerating power at the beginning of the interval time . M
The process of computation is now repeated to obtain Pa (n), (n+1) and (n+1). The solution in discrete form is thus carried out over the desired length of time normally 0.05 sec. Greater accuracy of solution can be achieved by reducing the time duration of interval.
Any switching event such as occurrence of a fault or clearing of the fault causes discontinuity in the accelerating power P a. If such a discontinuity occurs at the beginning of an interval then the average of the values of P a just before and just after the discontinuity must be used. Thus in computing the increment of angle occurring during the first interval after a fault is applied at time t = 0 becomes
2 ( t) 2 1 1 + + ( t) 1 = 0 + (Pa 0 + Pa 0 ) = Pa 0 (Because Pa 0- = 0) M M 2 2
If the discontinuity occurs at the middle of an interval, no special procedure is needed. The correctness of this can be seen from Fig. 15.
n th interval n - 1 th interval
n th interval n - 1 th interval
Fig. 15
Example 4 A 20 MVA, 50 Hz generator delivers 18 MW over a double circuit line to an infinite bus. The generator has KE of 2.52 MJ / MVA at rated speed and its transient reactance is Xd = 0.35 p.u. Each transmission line has a reactance of 0.2 p.u. on a 20 MVA base. E = 1.1 p.u. and infinite bus voltage V = 1.0 p.u. A three phase fault occurs at the mid point of one of the transmission lines. Obtain the swing curve over a period of 0.5 sec. if the fault is sustained.
E = 1.1 p.u.
G 20 MVA 50 Hz Delivers 18 MW Xd = 0.35 p.u. H = 2.52 MJ/MVA
V = 1.0 p.u.
Recursive equations are (n) = (n 1) + (n) Pre fault: During fault: 0.2 0.1 0.35 + 1.1 0.1 + 1.0 0.35 + 0.1 0.1 + 0.2 where (n) = (n 1) + 8.929 Pa (n-1) Pe =
1.1 x 1.0 sin 2.44 sin 0.45
X = 0.45 p.u.;
1.25 1.0
Initial calculations: Before the occurrence of fault, there will not be acceleration i.e. Input power is equal to output power. Therefore Ps = 18 MW = 0.9 p.u. Initial power angle is given by 2.44 sin 0 = 0.9; Pa 0- = 0; Pa
average
Thus 0 = 21.64
Thus (0.05) = 21.64 + 2.57 = 24.210 Subsequent calculations are shown below.
t sec. 00
+
deg. 21.64 21.64 21.64 24.21 31.59 42.89 56.87 72.30 88.28 104.44 121.02 138.90 159.65
Pmax
Pe
Pa = 0.9-Pe 0
8.929 Pa
0.88
0.324
0.576 0.288 2.57 4.81 3.92 2.68 1.45 0.55 0.18 0.426 1.30 2.87 2.57 7.38 11.30 13.98 15.43 15.98 16.16 16.58 17.88 20.75
0 average 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
.05
.1
.15 .2
.25
.3
.35
.4
.45
.5
Fig. 16
Swing curve, rotor angle with respect to time, for sustained fault is plotted and shown in Fig. 16.
Example 5 For the previous example, calculate the CCA and hence find CCT. Solution
cos C C
Ps = 0.9;
PS ( m 0 ) Pm 2 cos 0 Pm 3 cos m Pm 3 Pm 2
Pm1 = 2.44; Pm2 = 0.88; Pm3 = 2.0
0.9 ( 2.6748 0.3778 ) 0.88 cos (0.3778) 2 cos (2.6748) 0.47915 2 0.88
Thus CCA, CC = 118.630 Referring to the graph shown in Fig. 16, corresponding to CCA of 118.63 0, CCT can be obtained as 0.38 sec
SOLUTION OF SWING EQUATION BY MODIFIED EULERS METHOD Modified Eulers method is simple and efficient method of solving differential equations (DE) Let us first consider solution of first order differential equation. Later we shall extend it for solving a set of first order DE. The swing equation is a second order DE which can be written as two first order DE and solution can be obtained using Modified Eulers method. Let the given first order DE be
dx f ( t, x ) dt
(34)
where t is the independent variable and x is the dependent variable. Let (t 0, x 0) be the initial solution and t is the increment in t. Then t 1 = t 0 + t; t 2 = t 1 + t; t n = t n 1 + t
(35)
Thus (t 1, x 1(0)) is the first estimated point (t 1, x 1). Second and the final estimate of x 1 is calculated as x1 = x0 +
1 dx dx (0) |0 + |1 ) t ( 2 dt dt
(36)
where
(t 1, x 1) is now known. Same procedure can be followed to get (t 2, x 2) and it can be repeated to obtain points (t 3, x 3), (t 4, x 4) .. Knowing (t n - 1, x n - 1 ), next point (t n, x n) can be computed as follows: x n(0) = x n 1 + t n = t n 1 + t Compute
dx (0) dx |n which is computed at (t n, x n(0)). dt dt 1 dx dx (0) |n - 1 + |n ) t ( 2 dt dt dx |n 1 t dt
Then
xn = xn-1 +
(40)
Same procedure can be extended to solve a set of two first order DE. Let the DE be
dx = f1 (t, x, y) dt
and
dy = f2 (t, x, y) dt
y n(0) = y n 1 + t n = t n 1 + t Compute
Then x n = x n - 1 +
1 dx dx (0) ( |n - 1 + |n ) t 2 dt dt 1 dy dy (0) ( |n - 1 + |n ) t 2 dt dt
and
yn = yn-1 +
When per unit values are used and the machines rating is taken as base M=
H f
d = s dt d = K ( Ps Pe ) dt
d d generally of the form = f1 ( t, , ). However, now it a function dt dt d d of alone. Similarly, generally of the form = f2 ( t, , ). However, now it dt dt a function of alone.
Note that
Just prior to the occurrence of the disturbance, P s Pe = 0 and = s. The rotor angle can be computed as (0) and the corresponding angular velocity is (0). Thus the initial point is ( 0, (0), (0)). As soon as disturbance occurs, electric network changes and the expression for electric power Pe in terms of rotor angle can be obtained. During fault condition, Pe shall be computed by the said expression. Using Modified Eulers method 1 and 1 can be computed. Thus we get the next solution point as ( t1, 1, 1). The procedure can be repeated to get subsequent solution points until next change in electric network, such as removal of faulted line occurs. As soon as electric network changes, corresponding expression for electric power need to be obtained and used in subsequent calculation. The whole procedure can be carried out until t reaches the time upto which transient stability analysis is required.
Example 6 An alternator rated for 100 MVA supplies 100 MW to an infinite bus through a line of reactance 0.08 p.u. on 100 MVA base. The machine has a transient reactance of 0.2 p.u. and its inertia constant is 4.0 p.u. on 100 MVA base. Taking the infinite bus voltage as reference, current supplied by the alternator is ( 1.0 j 0.6375 ) p.u. Calculate the torque angle and speed of the alternator for a period of 0.14 sec. when there is a three phase fault at the machine terminals and the fault is cleared in 0.1 sec. Use Modified Eulers method with a time increment of 0.02 sec. Solution j 0.2 p.u. 100 MVA Ps = 100 MW H=4 E j 0.08 p.u.
I = (1.0 j 0.6375)p.u.
V = 1.0 0 0 p.u.
E = ( 1.0 + j0) + j 0.28 ( 1.0 j 0.6375) = 1.1785 + j 0.28 = 1.2113 13.3651 0 p.u. Initial rotor angle = 13.36510 = 0.2333 rad.
Initial point is: (0) = 0.2333 rad. (0) = 314.1593 rad. / sec.
= 0.2333 + ( 0 x 0.02 ) = 0.2333 rad. = 314.1593 + ( 39.2699 x 0.02 ) = 314.9447 rad. / sec. Second estimate:
d = 314.9447 314.1593 = 0.7854 dt d = 39.2699 ( 1 0 ) = 39.2699 dt
(0.02) = 314.1593 +
Latest point is: (0.02) = 0.24115 rad. (0.02) = 314.9447 rad. / sec.
= 0.24115 + ( 0.7854 x 0.02 ) = 0.2569 rad. = 314.9447 + ( 39.2699 x 0.02 ) = 315.7301 rad. / sec. Second estimate:
d = 315.7307 314.1593 = 1.5708 dt d = 39.2699 ( 1 0 ) = 39.2699 dt
(0.04) = 314.9447 +
Calculations can be repeated until the fault is cleared i.e. t = 0.1 -. The results are tabulated. Thus (0.1) = 0.4297 rad.; (0.1) = 318.0869 rad. / sec.
Once the fault is cleared, reactance between internal voltage and the infinite bus is 0.28 and thus generator out put is; Pe =
1.2113x 1.0 sin 4.3261 sin 0.28
In the subsequent calculation Pe must be obtained from the above equation. To calculate (0.12) and (0.12) Latest point is: (0.1) = 0.4297 rad. (0.1) = 318.0869 rad. / sec.
Latest point is: (0.1) = 0.4297 rad. (0.1) = 318.0869 rad. / sec.
= 0.4297 + ( 3.9276 x 0.02 ) = 0.50825 rad. = 318.0869 + (- 31.5041 x 0.02 ) = 317.4568 rad. / sec. Second estimate:
d = 317.4568 314.1593 = 3.2975 dt d = 39.2699 ( 1 4.3261 sin 0.50825 rad. ) = - 43.4047 dt
Thus
(0.12) = 0.4297 +
1 ( 3.9276 + 3.2975 ) x 0.02 = 0.50195 rad. 2 1 ( - 31.5041 - 43.4047 ) x 0.02 = 317.3378 rad. / sec. 2
(0.12) = 318.0869 +
First Estimate
t sec. rad. rad/sec d/dt d/dt rad. rad/sec d/dt
Second Estimate
d/dt rad. rad/sec
0.2333 0.2333 0.2412 0.2647 0.304 0.359 0.4297 0.4297 0.502 0.557
314.1593 314.1593 314.9447 315.7301 316.5155 317.3009 318.0869 318.0869 317.3378 316.3955 3.9276 3.1785 -31.504 -42.468 0.5083 0.5655 317.4568 316.4884 3.2975 2.3291 -43.405 -51.761 0.502 0.557 317.3378 316.3955 0 0.7854 1.5708 2.3562 3.1416 39.2699 39.2699 39.2699 39.2699 39.2699 0.2333 0.2569 0.2961 0.3511 0.4218 314.9447 315.7301 316.5155 317.3009 318.0863 0.7854 1.5708 2.3562 3.1416 3.927 39.2699 39.2699 39.2699 39.2699 39.2699 0.2412 0.2647 0.304 0.359 0.4297 314.9447 315.7301 316.5155 317.3009 318.0869
0 0.2333 13.37
SOLUTION OF SWING EQUATION BY RUNGE KUTTA METHOD Fourth order Runge Kutta (RK) method is one of the most commonly used methods of solving differential equation. Consider the first order DE
dx f ( t, x ) dt
Let (t m, x m) be the initial point and h be the increment in time. Then tm+1 = tm + h Fourth order RK method can be defined by the following five equations. xm+1 = xm +
1 ( k1 + 2 k2 + 2 k3 + k4) 6
where
k1 = f ( t m, x m ) h k2 = f ( t m +
k h ,xm+ 1)h 2 2
k h ,xm+ 2 )h 2 2
k3 = f ( t m +
k4 = f ( t m + h, x m + k3 )
Note that in this method, the function has to be evaluated four times in each step. Same procedure can be extended to solve a set of first order DE such as
dx f1 ( t , x , y ) dt
and
dy f2 ( t , x , y ) dt
1 ( k1 + 2 k2 + 2 k3 + k4) 6
ym+1 = ym +
1 ( 1 + 2 2 + 2 3 + 4) 6
1 = f2 ( t m, x m, ym ) h 2 = f2 ( t m +
k h , xm + 1 , ym + 1 ) h 2 2 2
k h , xm + 2 , ym + 2 ) h 2 2 2
k3 = f 1 ( t m +
3 = f2 ( t m +
k4 = f1 ( t m + h, xm + k3, ym + 3 ) h
4 = f2 ( t m + h, xm + k3, ym + 3 ) h
The initial solution point is ( 0, (0), (0)). When 4th order RK method is used, k1, 1, k2, 2, k3, 3, k4, 4 are computed and then the next solution point is obtained as ( t1, 1, 1). This procedure can be repeated to get subsequent solution points.
Example 7 Consider the problem given in previous example and solve it using 4 th order RK method. Solution As seen in the previous example, two first order DEs are
d = 314.1593 dt d = 39.2699 ( 1 Pe ) dt
Initial point is: (0) = 0.2333 rad. (0) = 314.1593 rad. / sec.
During the first switching interval t = 0+ to 0.1 sec. electric output power Pe = 0.
To calculate (0.02) and (0.02) k1 = (314.1593 314.1593) x 0.02 = 0 (0) + k1 / 2 = 0.2333; 1 = 39.2699 ( 1 0 ) x 0.02 = 0.7854
k2 = (314.552 314.1593) x 0.02 = 0.007854 2 = 39.2699 ( 1 0 ) x 0.02 = 0.7854 (0) + k2 / 2 = 0.2372; (0) + 2 / 2 = 314.1593 + 0.3927 = 314.552 3 = 39.2699 ( 1 0 ) x 0.02 = 0.7854
(0.02) = 314.1593 +
It is to be noted that up to 0.1 sec., since Pe remains at zero, constants 1 = 2 = 3 = 4 = 0.7854 To calculate (0.04) and (0.04) k1 = (314.9447 314.1593) x 0.02 = 0.01571 (0.02) + 1 / 2 = 314.9447 + 0.3927 = 315.3374 k2 = (315.3374 314.1593) x 0.02 = 0.02356 (0.02) + 2 / 2 = 314.9447 + 0.3927 = 315.3374 k3 = (315.3374 314.1593) x 0.02 = 0.02356 (0.02) + 3 = 314.9947 + 0.7854 = 315.7301 k4 = (315.7301 314.1593) x 0.02 = 0.03142 (0.04) = 0.24115 +
1 [ 0.01571 + 2 (0.02356) + 2 (0.02356) + 0.03142 ] = 0.2647 rad. 6
Latest point is: (0.02) = 0.24115 rad. (0.02) = 314.9447 rad. / sec.
To calculate (0.06) and (0.06) k1 = (315.7301 314.1593) x 0.02 = 0.03142 (0.04) + 1 / 2 = 315.7301 + 0.3927 = 316.1228 k2 = (316.1228 314.1593) x 0.02 = 0.03927 (0.04) + 2 / 2 = 315.7301 + 0.3927 = 316.1228 k3 = (316.1228 314.1593) x 0.02 = 0.03917 (0.04) + 3 = 315.7301 + 0.7854 = 316.5155 k4 = (316.5155 314.1593) x 0.02 = 0.04712 (0.06) = 0.2647 +
1 [ 0.03142 + 2 (0.03927) + 2 (0.03927) + 0.04712 ] = 0.304 rad. 6
Latest point is: (0.04) = 0.2647 rad. (0.04) = 315.7301 rad. / sec.
To calculate (0.08) and (0.08) k1 = (316.5155 314.1593) x 0.02 = 0.04712 (0.06) + 1 / 2 = 316.5155 + 0.3927 = 316.9082 k2 = (316.9082 314.1593) x 0.02 = 0.05498 (0.06) + 2 / 2 = 316.5155 + 0.3927 = 316.9082 k3 = (316.9082 314.1593) x 0.02 = 0.05498 (0.06) + 3 = 316.5155 + 0.7854 = 317.3009 k4 = (317.3009 314.1593) x 0.02 = 0.06283 (0.08) = 0.2647 +
1 [ 0.04712 + 2 (0.05498) + 2 (0.05498) + 0.06283 ] = 0.359 rad. 6
Latest point is: (0.06) = 0.304 rad. (0.06) = 316.5155 rad. / sec.
To calculate (0.1) and (0.1) k1 = (317.3009 314.1593) x 0.02 = 0.06283 (0.08) + 1 / 2 = 317.3009 + 0.3927 = 317.6936 k2 = (317.6936 314.1593) x 0.02 = 0.07069 (0.08) + 2 / 2 = 317.3009 + 0.3927 = 317.6936 k3 = (317.6936 314.1593) x 0.02 = 0.07069 (0.08) + 3 = 317.3009 + 0.7854 = 318.0863 k4 = (318.0863 314.1593) x 0.02 = 0.07854 (0.1) = 0.359 +
1 [ 0.06283 + 2 (0.07069) + 2 (0.07069) + 0.07854 ] = 0.4297 rad. 6
Latest point is: (0.08) = 0.359 rad. (0.08) = 317.3009 rad. / sec.
(0.1) = 317.1593 + 0.7854 = 318.0863 rad. / sec. At t = 0.1 sec., the fault is cleared. As seen in the previous example, for t 0.1 sec., electric power output of the alternator is given by Pe = 4.3261 sin
To calculate (0.12) and (0.12) k1 = (318.0863 314.1593) x 0.02 = 0.07854 Pe = 4.3261 sin (0.4297 rad.) = 1.8022 1 = 39.2699 ( 1 1.8022 ) x 0.02 = - 0.63 (0.1) + k1 / 2 = 0.4690; (0.1) + 1 / 2 = 318.0863 - 0.315 = 317.7713 Latest point is: (0.1) = 0.4297 rad. (0.1) = 318.0863 rad. / sec.
k2 = (317.7713 314.1593) x 0.02 = 0.07224 Pe = 4.3261 sin (0.469 rad.) = 1.9554 2 = 39.2699 ( 1 1.9554 ) x 0.02 = - 0.7504 (0.1) + k2 / 2 = 0.4658; (0.1) + 2 / 2 = 318.0863 + 0.3752 = 317.7111
k3 = (317.7111 314.1593) x 0.02 = 0.07104 Pe = 4.3261 sin (0.4658 rad.) = 1.9430 3 = 39.2699 ( 1 1.9430 ) x 0.02 = - 0.7406 (0.1) + k3 = 0.5007; (0.1) + 3 = 318.0863 - 0.7406 = 317.3457
k4 = (317.3457 314.1593) x 0.02 = 0.06373 Pe = 4.3261 sin (0.5007 rad.) = 2.0767 4 = 39.2699 ( 1 2.0767 ) x 0.02 = - 0.8456 (0.12) = 0.4297 +
1 [ 0.07854 + 2 (0.07224) + 2 (0.07104) + 0.06373 ] = 0.5012 rad. 6
To calculate (0.14) and (0.14) k1 = (317.3434 314.1593) x 0.02 = 0.06368 Pe = 4.3261 sin (0.5012 rad.) = 2.0786 1 = 39.2699 ( 1 2.0786 ) x 0.02 = - 0.8471 (0.12) + k1 / 2 = 0.5330; (0.12) + 1 / 2 = 317.3434 - 0.42355 = 316.91985 Latest point is: (0.12) = 0.5012 rad. (0.12) = 317.3434 rad. / sec.
k2 = (316.91985 314.1593) x 0.02 = 0.05521 Pe = 4.3261 sin (0.533 rad.) = 2.1982 2 = 39.2699 ( 1 2.1982 ) x 0.02 = - 0.9411 (0.12) + k2 / 2 = 0.5288; (0.12) + 2 / 2 = 317.3434 - 0.47055 = 316.87285
3 = 39.2699 ( 1 2.1825 ) x 0.02 = - 0.9287 (0.12) + k3 = 0.5555; (0.12) + 3 = 317.3434 - 0.9287 = 316.4147
k4 = (316.4147 314.1593) x 0.02 = 0.04511 Pe = 4.3261 sin (0.5555 rad.) = 2.2814 4 = 39.2699 ( 1 2.2814 ) x 0.02 = - 1.0064 (0.14) = 0.5012 +
1 [ 0.06368 + 2 (0.05521) + 2 (0.05427) + 0.04511 ] = 0.5558 rad. 6