.A DIGITAL COMPUTER SOLUTION OF SYNCHRONOUS MACHINE DIFFERENTIAL EQUATIONS IN SOME POWER SYSTEM STABILITY PROBLEMS By Ralph Wayne Gilchrist AN ABSTRACT submitted to the School of Advanced Graduate Studies of Michigan State University of Agriculture and Applied Science in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical Engineering 1960 Approved 4 1- 4w if". .fitfi ABSTRACT It is common practice to make simplifying assumptions in order that convenient or familiar means may be used in the study of synchronous ma- chines. For the purpose of this thesis, these assumptions will be con- sidered in two groups. With the first group of assumptions, there re- sults the familiar system of non-linear differential equations interre- lating voltages, currents, speed, and torque for the synchronous machine. This type of system of equations is not generally amenable to solution except by numerical methods. Rather than attempt to solve this formid- able system of equations directly, engineers working in the area of power system stability use a second group of simplifying assumptions in order to solve the system of equations by convenient means. Using both groups of assumptions, the conventional method thus developed over the years is to solve part gf the equations by modified steady state techniques and then apply numerical methods to the remaining equations in order to obtain a solution. Similar techniques are used by engineers in other areas where the transient operation of synchronous machines is of interest. The study of a synchronous motor with a periodically varying load or with impact loading is an important example of these areas. Recently, much work has been done toward applying‘the digital comp puter to the study of the power system stability problem. However, the trend is toward programming the conventional techniques with their two sets of simplifying assumptions and empirical methods. In contrast, this thesis is a report on an investigation of the numerical solution of the complete system.of synchronous machine non-linear differential equations without the second group of simplifying assumptions. -2- Two classes of power system stability problems are studied here and it is found that for the cases considered: 1. The solution can be carried out in detail without the second group of simplifying assumptions. 2. The solution can be carried out under transient conditions with- out developing or using extensive empirical relations that are necessary in using modified steady state techniques for the transient prOblems. 3. The form of the results reported here gives more stability in- formation and gives the information more directly. h. The direct solution of the complete set of synchronous machine- equations rather than solution by conventional means, produces results that add to and verify accepted theory in some cases and contradict ac- cepted theory in other cases. A conclusion from the investigation reported here is that, for the class of problems considered in this thesis, methods developed before the advent of the digital computer should not be used as the basis for digi- tal computer investigation. The more fundamental synchronous machine equations can now be solved directly. The physical and mathematical structure for the cases considered here is the same as the structure in other areas involving transients in synchronous machines. Therefore the advantages of direct solution of the synchronous machine equations extend to areas other than power system stability. A DIGITAL COMPUTER SOLUTION OF SYNCHRONOUS MACHINE DIFFERENTIAL EQUATIONS IN SOME POWER SYSTEM STABILITY PROBLEMS By Ralph Wayne Gilchrist A THESIS Submitted to the School for Advanced Graduate Studies of Michigan State University of Agriculture and Applied Science in partial fulfillment of the requirements for the degree of DOCI‘CB OF PHILOSOPHY Department of Electrical Engineering 1960 7._‘_..» ACKNOHLEDGMENTS The author wishes to thank Dr. M. B. Reed, his thesis advisor, and Dr. H. E. Kbenig, his maJor professor, for their assistance and encouragement during the preparation of this thesis. TABLE OF CONTENTS ACICNMAEDWTS I I I I I I I I I I I I I I I I l I I I I I I I CHAPTERS I. II. III. IV. “I INTRODUCE I ON I I I I I I I I I I I I I I I I I I I I I I THE ROUND ROTOR SYNCHRONOUS MACHINE EQUATIONS AND THE TRANSFORMATIONS OF THE VARIABLES. . . . . . . . . . . . . TWO POWER SYSTEM STABILITY PROBLEMS . . . . . . . . . . . 3.1 The First Class of Problems . . . . . . . . . . . . 3.2 Some Refinements on the Rotor-Stator Coupling Terms 3.3 The Second Class of Problems . . . . . . . . . . . . SOLUTION OF A STABILITY PROBLEM FOR A SYNCHRONOUS MACHINE ON AN INFINITE BUS . . . . . . . . . . . . . . . h.l The Steady State Initial Condition Calculations . . h.2 Numerical Study of the Stability Problem of a Syn- chronous Machine on an "Infinite Bus" . . . h.3 Standard Procedure for Solving the Stability Problem of a Synchronous Machine on an "Infinite Bus" . . . h.h Comparison of Results and Methods . . . . . . . . . Two MACHINE SYSTEMS I I I I I I I I I I I I I I I I I I I 5.1 The Differential Equations of the Two Machine System 5.2 Steady State Initial Conditions for the Two Machine Problem . . . . . 5.3 Numerical Study of the Stability Problem of a Two Machine System . . . . . . . . . . . . . . . . . . . 5.h Standard Solution for the Two Machine Problem . . . 5.5 Comparison of Results and Methods . . . . . . . . . WY. I I I I I I I I I I I I I I I I I I I I I I I I APPENDIX A APPENDIX B BIBLIOGRAPHY (1i) Page (1) BB 24 28 39 A3 1&5 51 67 I. INTRODUCTION The purpose of a study of a synchronous machine is either to obtain design information or to predict operating conditions for a machine. In either case, use is made of the relation between the machine parameters: inductances, resistances, inertia and the machine variables: voltages, currents, speed, torque. The design problem is the problem.of obtaining the proper machine parameters for a desired range of the variables. The problem of predicting operating conditions is the problem.of obtaining a subset of the variables for conditions imposed by specifying the machine parameters and the remaining variables. To aid in the study of a synchronous machine, it is desirable to ob- tain a mathematical interrelation of the machine variables. Laboratory meter indications show the relationship among the variables to be so com- plicated that either exact mathematical interrelations cannot be found or the resulting mathematical system.of equations cannot be solved. Some simplifying assumptions, which lead to a mathematical system.of equations that can be solved, have become conventional. The system of equations resulting from the simplifying assumptions usually gives results that have acceptable correlation with laboratory meter indications over a limp ited range of the variables and for limited waveforms of the variables. One of the prdblems encountered in Obtaining a mathematical relation for the variables, is the irregular waveforms of the instrument indica- tions. This prdblem is commonly treated by considering only the lower frequency harmonics. Often only the fundamental frequency terms are used in the representations. Another prdblem is encountered in the fact that the apparent induc- -2- tance and resistance coefficients are not constants but are rather com- plicated functions of currents, power factor, field saturation, and other factors. There are clearly defined procedures for determining inductance and resistance numbers which are mean constants that partially take into account the dependence of the coefficients on current, power factor, and saturation. The fundamental harmonic representation and constant inductance and resistance coefficients are taken here as the first of two groups of sim— plifying assumptions. With this first group of assumptions, there results the standard system of differential equations for synchronous machines non-linear in speed torque and current variables. In order to solve the synchronous machine equations by convenient or familiar means, it is con- ventional to make a second group of simplifying assumptions also. The second group of assumptions are: constant speed operation, steady state sinusoidal armature voltages and currents, constant direct current field. With both groups of assumptions, the equations reduce to an algebraic system of equations which are readily solved. Relaxing any of these restrictions would yield pertinent information and give mathematical results which more closely correspond to laboratory observations of the meter indications. Much work has been done in the study of the effects of higher harmonics and in the study of the effects of saturation, current, and power factor on the inductance and resistance coefficients. Elaborate techniques have been developed to make steady state, constant speed techniques yield results correlating with observa- tions under transient conditions. These techniques usually involve sep- arating the equations into electrical and mechanical sets. The electrical -3- equations are solved by modified steady state methods and the mechanical equations are solved as differential equations, usually by numerical methods. Little has been done, however, in the area of solving the com- plete set as simultaneous non-linear differential equations. The pre- sently accepted standard techniques were developed before the advent of the digital computer, and numerical methods for solving merely the me— chanical equation as a non-linear differential equation involve long and tedious calculations. With the modern high speed computers, however, there appeared, at the beginning of the research reported.here, to be no reason why the complete set should not be solved as non-linear differen- tial equations by numerical methods. This thesis demonstrates the technique for solving the complete set of equations by numerical methods by applying the techniques to two classes of power system stability problems. Not only are the prdblems solved without making the assumptions of constant speed, constant field current, and sinusoidal armature voltages and currents, but for the classes of prdblems considered the solution process is simpler than the standard procedure. The form of the answers gives more stability infor- mation and gives the information more directly. Further, the greater generality of the techniques used in this thesis yields information which in some cases adds to or verifies the accepted theory and in other cases contradicts accepted theory. In particular, this additional information was obtained relative to the variation of the field current and the vari- ation of the direct and quadrature components of armature currents under transient conditions of operation. The techniques used here in solving the systems of equations for the -u- two classes of power system stability problems can be considered in three parts: First, the variables are changed in order to obtain a convenient form for the coefficients. The variables used are closely related to the variables of A. Blondel'sl two reaction analysis or the variables resulting from a transformation of variables illustrated by R. H. Park.2 Secondly, the differential equations are manipulated into standard form to permit numerical solution. In order to assure a unique continuous solution, the form of the equations must be such that they can be put into the standard form. Thirdly, an existing digital computer library routine is used to solve the system of equations. The routine used is based on the fourth order Runge-Kutta formulas. In conjunction with the differential equation solving routine, a master and a function routine must be supplied. The master and the function routines depend on the particular problem to be solved and on the amount of data the computer is required to print out. The greater simplicity of the technique of this thesis is evident, both in theory and application. The presently accepted standard tech- niques for solving the machine equations use the simplifying assumptions of constant speed and sinusoidal armature voltage and current variables in part of the equations. The constant speed assumption makes it possible to consider the equations in two separate groups. One group, called the electrical equations, is a system of linear differential equations con- taining only electrical variables. The other group, called the mechani- cal equations, or the swing equations, is made up of one equation for each machine under consideration. The swing equations are differential equations non-linear in speed, field current, and armature current. The -5- further assumption of sinusoidal steady state armature voltages and cur- rents and a change of armature variables to the familiar direct and quadrature components reduce the electrical equations to an algebraic system of equations. Under steady state conditions, if the power, power factor, speed, and armature voltages are specified, the algebraic system of electrical equations can be used to determine armature and field cur- rents. Neglecting losses, the armature and field currents determined from the electrical equations satisfy the mechanical equations. But, under fault conditions, the electrical and mechanical variables are in a transient state and the simplified equations cannot be expected to yield a solution which correlates with a solution of the differential equations. To compensate for the difference in differential equation solution and a steady state solution, the standard solution proceeds as follows: 1. Initial armature currents, field currents, and angular positions are determined from the steady state electrical equations. 2. New reactances, called transient reactances, are used with the results of the steady state solution, previously Obtained, to determine voltages and voltage angular positions. These voltages are called the voltages behind the transient reactances. 3. It is assumed that the voltages behind the transient reactances and the transient reactances of the machines do not change with the changes in speed and current to be considered. A. The voltages behind the transient reactances and the transient reactances do change as fault and the fault removal change the network parameters, however. -6- 5. Using initial conditions defined by previous steady state calcu- lations and by the type of fault, the swing equations are solved.by common numerical methods. For each increment of time assigned in the numerical process, corresponding increments of changes in phase angles for the voltages behind the transient reactances are determined. These changes in angles are incorporated into the electrical power terms in the swing equations when the next time increment is applied to the solu- tion of the swing equations. 6. The solution is continued by adding increments of time until sufficient information is available on the phase angles of the voltages behind the transient reactances to determine if the system is stable* or unstable. 7. Since the electrical power is calculated in terms of transient reactance, voltage behind the transient reactance, and the phase angle of the voltage behind the transient reactance, different power formulas must be used for different states: steady state, fault on, fault removed. 8. If the changes in field currents are to be considereda, the voltages behind the transient reactances must be modified for each in- crement of time considered. The synchronous machine equations are solved directly in this thesis with only the first group of simplifying assumptions. Thus, the solution under transient conditions is obtained without formulation in terms of the so-called "transient reactance" and "field variation" effects and so iA system is considered unstable if any of the synchronous machines fall out of step with the system. -7- ‘without defining processes in which these reactances are used to give results correlating with observations under transient conditions. A common power system stability prdblem is that of determining how fast the circuit breaker must operate and remove a faulted line for the system to remain stable. Conventionally, the critical interval is ob- tained by a criterion which determines a critical rotor displacement as the variable. In order to determine the critical time, the numerh:al solution, by standard methods, must proceed up to the critical rotor dis- placement angle, although it need not be carried out until instability is indicated. An alternative method, illustrated by Kimbarka, uses 232? calculated swing curves. The use of pre-calculated swing curves requires the steady state, constant speed, constant field current assumptions. The pre-calculated swing curves cannot be used for some limiting cases, and the pre-calculated swing curves do not represent conditions after the fault has been removed. The techniques of this thesis do not require the assumptions of either the criterion or the pre-calculated swing curve techniques. Further, the variables determined by the general solution presented here include both time and angle variables for the fault inter- val and for the interval after the fault has been removed. The switching time and the type of fault are controlled by the master routine in the computer program. The results reported here of a typical problem indicate that the field current variation, commonly given considerable attentiong, is in- significantly small even under extreme conditions, whereas the commonly Jaeglected change in speed is really a more dominant variation. The re- :sults of the numerical solution for the differential equations indicate ‘— -8- that the direct and quadrature currents of Blondel's analysis or Park's transformation have characteristics much different under transient condi- tions than that usually assumed in standard steady state analysis and that usually assumed in modified steady state analysis applied to the power system.stability prOblem. II. THE ROUND ROTOR SYNCHRONOUS MACHINE EQUATIONS .AND THE TRANSFORMATIONS OF'THE VARIABLES In this chapter the synchronous machine equations, obtained with only the first group of simplifying assumptions, are listed and all co- efficients are defined. Next, the transformations are listed and the equations in the new variables are given. The round rotor synchronous machine variables that are Observed by means of meter indications are: the three-phase voltages of the machine armature, the three-phase currents of the machine armature, the field voltage and current, and the speed and torque at the shaft. With the usual simplifying assumptions that the interrelationship of the variables is represented with sufficient accuracy by considering only the funda- mental harmonics and by considering the resistances and the inductances as constants, the mathematical interrelations of the variables are: Wt) 26’ +1.2” ggat’srm 480.) (2.1)a = as dt 65 vf(t) §E.Z’m(¢) IRE-5331'“ if“) _ l 9 83 8" 3 + (2-1» T--2 Us“) 114”] 75 arms) L.» if“) (B + J-gr 5 dt where (2/3“) and 005”) are column matrices representing the line to -9- neutral voltages and currents of the three stator (armature) phases. Specifically, —"a(*'fl 380$ 7/8“) vb(t) <93“) = ib(t) vcm , 1cm The functions vf(t) and if(t) represent the rotor (field) voltage and current respectively. The functions T(t) and.¢(t) are the shaft torque and position respectively. The coefficient matrices fless and ‘Qfss are of the form 72A 0 9-1 _RA o 0—7 WS 3 = 0 BB 0 = o R A o L_o o 129— Lo 0 R1}. ,. _ LAA LEA LG: LAA LAB LA; ogss = LAB LBB L013 = LAB LAA LAB :AC LBC I‘cc_J L __LA.'B LAB LAA_ where LAA, LBB, and LCC are the self inductance constants for the respec- tive stator phases and LAB’ LBA’ L , LC , LBC’ and LCB are the mutual inductance constants representing the coupling between stator phases. The constants Rff and Lff represent the resistance and inductance coefficients of the rotor. The coefficient matrix, caf;r(¢), is of the form HI cos 9 er cos (a - 1209) 0881190 = align) = sr r cos (9 - 2h00) let-- a: ‘—v .. 4'11‘1‘ It‘lfllll \ i: -10- = 2 with e 2 ¢ and p representing the number of poles of the machine and 6 representing the displacement of the rotor in electrical angle measure. The coefficients B and J represent the mechanical damping constant and polar moment of inertia respectively. A change of stator variables greatly simplifies the interrelation- ship of (2.1). variables closely related to those to be used here were develOped by A. Blondel1 in order to more accurately consider saliency effects of salient pole synchronous machines. R. H. Park2 illustrated that the new variables can be considered as a set of variables resulting from a mathematical transformation of variables. Park's transformations on the stator variables are: Hohfl _1 1 1 _ an)“ (2.2) vd(t) = ~13- 2 cos a 2 cos (av-120°) 2 cos (o+120°) vb(t) quh)‘ _:2 cos a -2 sin (ca-120°) -2 sin ($1203 ch(t)d _io(t)- _ 1 1 1 _ Hath.)-1 i am = § 2 cos a 2 cos (9-1200) 2 cos (9.120% in“) Liq(1-.) -2 sin e -2 sin (9-120°) -2 sin (9+120"_) bic(t)_ It is a formidable task to apply the transformations of (2.2) to the equations of (2.1). This application of the transformations would be necessary, however, in order to rigorously define the coefficients in the :resulting equations in terms of the observed relationships (2.1). Koenig3 has shown these transformations (2.2) to be the product of 'three non-singular transformations. The three component transformations applied to the stator voltage variables are: -11- v:(t)— l l l va(t)— (2.3) v;(t) 71; 1 e3a/3 633/3 vb(t) L";("')_ 4.1 ewes/3 east/3 30m- 2:“)... ‘1 0 o- 1:“; (2.1.) v:(t) = o {39 o v;(t) c:(t)4 LO 0 eJ€__ [Jr's-(t).-J where 92(t) is the conjugate of v:(t). -}3tf1 .2 o ow'lgsf (2.5) vd(t) 31—2 0 1 1 v:(t) Ab :9“)... _0 J -1_ Us“). The transformations defined in (2. 3), (2.4), and (2.5) utilize symmetries and inherent characteristics of the coefficient matrices and with theorems developed by Koenig and Blackwellh, these transformations are relatively easy to apply. The transformation defined in (2.3) is the familiar symmetrical component transformation. Application of the transformation of (2.3) followed by the transformation in (2.h) defines the backward sequence variables. Application of the transformation in (2.5) to the backward sequence variables separates the real parts of the backward sequence com- plex variables from the imaginary parts of the backward sequence complex variables. The new armature real variables resulting from the three transformations are essentially the same as the direct-axis and quadra- ture-axis components of armature voltages. Even though salient pole synchronous machines, for which these new 'variables were developed, are not considered in this investigation, the -12.. new variables are in a convenient form for numerical solution of the differential equations. After application of the transformations defined in (2.3), (2.h), and (2.5) for the voltage variables of (2.1) and the same transformations are applied to the current variables of (2.1), the machine equations be- come: 0 o o d o vs(t) RS+LS at o o o 15(t) + + d ° + d . (2 6)a Va“) _ ° Rs+Ls dt Ls JELsr dt 1d“) ' ‘ ° + + + d ' . vq(t) o - 9 Ls RS+LS at -ofensr iqa) d d Vf“) o JELsr 317. o Rf+Lf at if(t) f2 (2.6» T(t) = ? LsrP1q(’°) if Figure 3.1 I At any particular operating point, i f = I f, the voltage-current relation is given by the slope of the secant line through P. Thus, the éJELsr terms in (2.6)a and (3.1) can be assigned a value equal to the slope of this secant line. For incremental changes about I f, however, the Oper— ating point moves along the tangent line. This refinement is incorporated into the term v(t) = évfingr if(t) by defining if(t) = If + if1(t) where I is a constant or reference value of if(t) 16' 1f1 If .eJ—ELsr is the slope of the secant line in Fig. (3.1) and BK is is the time varying component of if(t) the slope of the tangent line at if(t) = If, then v(t) = 9~f2Lsr If + ex lfl(t) .All of the terms in (2.6)a and (3.1) which involve field and armature -19- coupling are appropriately modified. The machine equations become 7.4a _ 32%: ads 5L? K at — 3.05)" (3.2)a vq(t) + ism/2LBr If = - ' L; 324i; dg't' - ax iq(t) vf - Rf If K 3‘11; 0 Rf+Lf g; in“) __ .L 1.. .1 c _ . . d2 (3.2)b Pin = e.~f2Lsr If iq(t) + ex ifl(t) iq(t) + M d? e The value of the slope of the tangent line changes with a change in operating point and a given slope would be in error for a large increment of change in current. The use of the tangent slope is an improvement over the use of the secant slope if if1(t) remains small. A typical value of slope OK for a synchronous machine at a typical operating point is about 0.6 or 0.7 p.u. with b = 377. This value Of OK as compared with typical values of 1.2 to 1.8 for b~f§isr shows the sig- nificance of the refinement on the field decrement effect. All terms in- volving field.and armature coupling are appropriately modified before the problem is programmed. 3.3 THE SECOND CLASS OF PROBLEMS The second class of prOblem considered here is the stability study for two interconnected machines. The significance of this class of prOb- lem lies in the fact that this prOblem represents an extension of the technique of the first prOblem toward a multi-machine system, and in the fact that with certain assumptions two interconnected systems fall into this class. Examples of problems of this class include: 1. Two synchronous machines are connected by a transmission line. One machine is operating as a generator and the other as a -20- motor. For a sudden increase in load on the motor, determine whether the machines will fall out of step or not. 2. Two power systems, each represented by an equivalent syn- chronous machine, are interconnected by a transmission line. For a fault on the transmission line which is cleared in a specified time, determine whether the systems remain in synchronism or not. 3. Two power systems, each represented by an equivalent syn- chronous machine, are interconnected by a transmission line. Sys- tem.A supplies power over the line, and system B receives power from the line. For a sudden increase of demand in system B, or for loss of a generating unit due to local fault, system.B require- ments become P + AP. Determine whether the two systems will remain in sychronism or not. under the assumptions that permit a system to be represented by an equivalent synchronous machine, the three prOblems are mathematically identical. The first of the three is used in this thesis to illustrate the technique. PROBLEM: For the second problem of this thesis, consider a synchro- nous generator supplying 0.5 p.u. power over a transmission line to a synchronous motor receiving 0.5 p.u. power. For a sudden change in load on the motor from 0.5 p.u. to 0.6 p.u., will the machines stay in syn- chronism? Repeat for the case where the motor load changes from 0.5 p.u. to 0.8 p.u. power. Initially, let the voltage at the motor terminals be 1 p.u. volts and the motor power factor be 0.85 lag. The machine parameters are assigned typical per unit values as in the first prOblem. The direct axis synchronous reactance of the motor -21- is taken as 0.8 p.u. and the H constant of the motor is taken as 2. Let the motor phase resistance be 0.005 p.u. and let the b‘IELsr be 1 p.u. for both the motor and the generator. Assign the generator plus the transmission line the direct axis synchronous reactance of 1.2 p.u. per phase and a phase resistance of 0.005 p.u. Let the generator H constant be 6. IV. SOLUTION OF A STABILITY PROBLEM FOR A SYNCHRONOUS MACHINE ON AN INFINITE BUS This chapter presents the standard and the thesis techniques in general and for the particular problem. Results are presented and com- pared. Methods are compared and the extension of methods to other prob- lems and other areas is discussed. h.l THE STEADY STATE INITIAL CONDITION CALCULATIONS The pre-fault conditions are the constant speed, steady state, and constant field conditions under which equation (2.6)a reduces to vdw Rs wLB 0 1d (11.1) v = - 0113+ R+ - (in/EL i q s 5 er q va- —0 0 Rf J L1f_ The pre-fault conditions given in the statement of the first problem are: (VI = 1 p.u. voltage at the bus. P = 0.5 p.u. power at 0.85 power factor lag. wL: = 1.2 p.u. reactance R: = 0.005 p.u. resistance [II = .5882 p.u. current In addition to equation (T1), the following equations apply for -22- steady state sinusoidal conditions: vd(t) + [V] sin 5 vq(t) p Q: A ¢+ V II - |I| sin (5 + s) |Ii cos (6 + B) - |v| cos 5 i (t) q where B is the power factor angle 8 is the phase angle between the bus voltage, V, and the field excitation‘voltage, Ef = (in/ELM. If Usually, this steady state problem is solved neglecting resistance and the various quantities are pictured on a vector diagram as in Fig. (h.1). Figure h.l Solving (h.l), the results are: Ef = 1.5 /90 + 23.60 p.u. volts 5 = 23.6° Id = -0.h96 p.u. current I = 0.335 p.u. current -23- If = 1.5 p.u. since Ef = If in p.u. vd = 0.h0 p.u. volts vq = -0.916 p.u. volts Further, since the pre-fault conditions are considered as constant syn- chronous speed, the initial values of 9 and O are: 9 377 radians per second 6 wto + s + ./2 = s + 1/2 = 0.1+116 radians In steady state terms, the angle 8, sometimes called the torque angle, is the angle represented in Fig. (h.l) as the phase angle between the bus voltage and the voltage generated in the armature due to the field current. In equations (2.1), and in the subsequent transformations, 8 is the angle of deviation of the rotor from the no-load synchronous speed position. A synchronous machine is unstable if 6 increases with- out limit. If a system is stable, the variation of 5 is generally an oscillatory function of time. If a machine is stable for the first swing, it is classified as stable, for by the end of the first swing, regulator action and prime mover governor action will have begun and these actions tend to stabilize the machine. It is generally accurate to assume that prime mover governor action does not take place until after the first swing of the oscillatory 8 and does not enter into the stability calculations. Therefore, it is accurate to assume that the prime mover power input remains constant during the first swing. Some authors7’8 have developed formulas to take into con- sideration generator field control by regulator action. The regulator action is not usually initiated before the fault has been on for 0.2 or 0.3 seconds, however. ' —-- h-C‘r—v—‘n: \. -214... h.2 NUMERICAL STUDY OF THE STABILITY PROBLEM OF A SYNCHRONOUS MACHINE ON AN "INFINITE BUS" by: and The mathematical statement of the problem to be solved here is given (1) equations (3.2)a and (3.2)b, (2) the initial conditions defined in section n.1, (3) the fault conditions, (4) the switching conditions. The normal form of the system equations given in (3.2)a and (3.2)b as required for a numerical solution by the Runge-Kutta method is: (h.2) V'— d d agect) fish) L. 51d“) d . Etiq(t) d Rim“) — [— + e + 1 K2 -‘ [vd(t)'Rsid(t) ' aLsiq(t)]‘L%’- L;2(Lf ‘ Lg) [ [vd(t)-R;id(t) - 6L; iq(t)] [vq(t) - R; iq(t) - OJ-ZLBrIf + b L:id(t) + bx ifl(t)] ii: _ _ + _ - + K i _ R L [vd(t) R8 id(t) 9 L15 iq(t):] L31}: - SJ f 1fl(t)/ f b [pm - 5 (J2 Lsr If + x ifl(t)) iq(t)] i- where, from.the transformations defined by equations (2.3), (2.h), and (2.5) and vd(t) vq(t) 6(t) IV, sin 5 - [V[ cos 5 am.+ 5 + 1/2 The flow diagram for the computer program Which was used to solve these equations is shown in Fig. (h.2). -25- Stavcd StOYta Data. pat“ +\ H >>** f:'- (a /' Nosfev fir -—-6— -<- +A V A 1".5' 53h 13.1 <10E> ‘Priht ou.t Figure h.2 The digital computer library routine, F-6, was used to solve the dif- ferential equations. The F-6 routine was chosen because it uses the fourth order Runge-Kutta10 formulas. The truncation error for the fourth order Runge-Kutta formulas is of the order of (At)5, where At is the increment of time used for each step. .Also, the F-6 routine limits the amount of error by an automatic control of the size of the increment. Both of these, the fourth order Runge-Kutta formulas and the increment control, tend to lnake the program slow, but both are valuable assets when the exact char- acter of the variables being investigated is unknown. The differential equation solving routine, the master routine, and tflae function routine are described in greater detail in Appendix B. With the program described in Appendix B available, the solution to -26- the first problem is obtained simply. Merely supply the program with the specified initial conditions, the parameters for the particular machine, and master routine counters. The scaled initial values, scaled coefficients, and constants used in the problem.are listed. The numbers in memory location 3 through 8 as required by F-6. In 3 In A In 5 In 6 In 7 In 8 In 10 the integer 30 indicates the location of the first of the sequence of initial conditions. the integer 36 indicates the sum of entry 3 and n, the num- ber of equations in the system. the integer 13 indicates the scale factor on the fr, 2'13. the integer 10 indicates the time increment, 2-10 see. the integer 36. The yr will have their error held to about 1 in the 36 binary bit. the integer 160 gives the memory location for the entry into the fr routine. the integer 100 gives the memory location of the first of a sequence of scaled coefficients and constants. The initial values of all the variables are scaled by 2'13 and read in as follows: Location Variable Scaled value 30 t 0 31 .9 = b 0.01.602 32 e 0.0002h20 33 1f1(t) 0 3“ 1d(t) -o.oooooosu6 35 iq(t) o.ooooh0882s The coefficients and constants for the fr routine: Location 100 101 102 103 10k 105 106 107 108 109 110 111 112 113 11A 115 116 Quantity 1/ J—zLBr M K + /L,3 1 /L; 1 /Lf P/M V Work space Work space Work space Work space 8 “/2 5/,“ Scale 2'9 2'3 2-13 20 Scaled Value 0 .73633 0.01592 0.500 0.61360 0.1666 0.0038910 0 0.001h6lI-8 0.0015915 0 . 00021+h1h program 0.0001917h7 0.63662 For master routine, read the following counters to control read out and.switching. If the integer k is placed in the 29th word of the master routine, :read out will occur after each (k+1)st calculation. If the integer p is placed in the 36th word of the master, switching «ar fault removal occurs after p read outs. If the integer q is placed in the 38th word of the master routine, —28- the program completes q read outs before stopping. The master reads out t, b, 6, ifl(t), 8, id(t), iq(t). However, if zero is placed in the 33rd word of the master routine, only t, 5, and 0 will be read out. The pre-set parameter in location 3 through 8 must be read in before F-6 is read in. All the other information can be placed at the end of the program. At the end of the read in, transfer control to the left side of the first word of the master routine and the computer gives the results. The critical switching time can be determined by altering the counter which is the 36th word of the master routine. The solution was Obtained for various switching times converging on a critical value. Some of the results appear in Graphs n.1, n.2, h.3. h.3 STANDARD PROCEDURE FOR SOLVING THE STABILITY PROBLEM.0F A SYNCHRONOUS MACHINE ON AN "INFINITE BUS" If a numerical solution is obtained for the system equations as in section (he), the technique is the same for all single mchine stability problems. In sharp contrast, standard.methods use a variety of techniques for various types of problems and also utilize additional assumptions and formulations to obtain results. In the interest of completeness and con- trast, section (h.3) presents some of these standard techniques. The standard procedure for the solution of the stability prOblem of a synchronous machine on an infinite bus is started.with the steady state solution of section (h.1). The values and.phase positions of the initial q’ vd’ vq, Ef and 6 are determined. I Usually, the equation (3.1) is written in the form d2 (3.1) P = m 2L:5r if(t) iq(t) + M 5-2— 95 id, 1 hoaxes“. S N run/pm w _p 8.0 506 86 ottoman. S N W06 v0.0 MOO N66 3.0 Q6 7rd 01 fuzz/”3 . ha “x .«u\. \QMNUM. Qx MN. \ '31- -32- By neglecting the resistance, the electrical power term can be written ‘ II I V E _ . . _ f . (14.3) Pelect — mJ-2Lsr 1f(t) 1q(t) - ---——-xd 81D 5 where IV! is the p.u. bus voltage lEfI is the p.u. voltage component due to the stator direct current field 8 is the torque angle + X d = 0 LS the direct axis synchronous reactance Standard methods, in order to use equation (4.1) with transient con- ditions, contain some modifications for equation (h.l) coefficients and contain some procedures for using equation (h.l). For a given steady state armature voltage, power, and power factor, equation (h.l) can be used to find v d is modified by replacing1 IE! = [sq + ,3 Ed. is assumed to remain constant during the transient interval if the field decrement is neglected. During this interval if losses are neglected, equation (1+ . 3) becomes . mug: _ ( .5) de sin 8 - Pelect Since 6 = mt + n/2 + 6, equation (3.1) becomes v E 2 (11.6) P =L-J—L‘Llsinb+M-‘Lb in X dT dt2 Under the assumptions necessary to obtain equation (11.6), a solution to (1.6) will indicate stability or instability. Standard techniques for solving a stability problem use a numerical solution of (1+.O) to determine -3h- stability characteristics. For this particular prOblem.of a short circuit at the bus, coeffi- cients in equation (h.6) are obtained using steady state and modified steady state methods. Then, from the time of occurrence of the fault until the circuit breaker removes the fault, equation (h.6) is solved with V set equal to zero. From.the time of the circuit breaker switch- ing and continuing to the end of the calculations, V is returned to the infinite bus status. The solution is carried out until the value of the angle 5 indicates stability or instability. Various circuit breaker oper- ating times must be considered until the time is found for which the ma- chine operation becomes unstable. This time is called the critical switch- ing time and the corresponding angle is called the critical switching angle. For any switching time less than the critical value, the machine operation would be stable, and for any switching time greater than the critical value, the machine operation would be unstable. Some important variations on the direct solution of equation (h.6) are now considered. In order to avoid carrying out the numerical solution of equation (4.6) for various switching times until a critical switching time is deter- mined, a criterion has been developed for determining the critical switch- ing angle. Except in simple cases, the determination of this critical angle requires a trial and error graphical integration approach. Further, even with the critical angle known, the switching time must be determined by numerical solution of equation (h.6) out to the critical angle. For the particular problem of a short circuit at the bus as in this first the- sis problem can be solved most readily by the equal area criterion. The -35- solution is presented in the following paragraphs. Under the appropriate assumptions, equation (h.3) was obtained. V E _ f (h-3) Pelect - ———i;r- sin 5 A sketch of P as a function of 8 is given in Fig. (h.h). The elect horizontal line represents the initial power supplied to the generator, 80 the initial torque angle, 50 the switching angle. The angle, 5c, is the critical switching angle when the cross-hatched area above the initial power line is equal to the cross-hatched area below the initial power line. —;> 03% /% ///// l 25 62: 1+ Figure h.h For the initial steady state condition [VI = 1 p.u. Pr) = 1.5 p.u. x(1 = 1.2 Pelect = 0.5 p.u. so = 23.6° = .h116 radians P : gl)§1.§) elect 1.2 for 5c the critical value 180° - 23.60 5 1.25 sin b d 8 - .5 (1800 - 23.6O - 5c) = .5 (5c - 23.6°) sin 8 = 1.25 sin 5 c then -36- + 0.01h O 0 m 0? 0 II 0 0’ ll 39 Next, with IVI equal to zero for short circuit at the bus, the criti- cal 6c is determined to occur at 0.28 seconds by numerical solution of equation (4.6). .A second variation from the direct solution of equation (h.6) becomes necessary if a study is to be made of the field decrement effect and/or if a study is to be made of the incrementing of the field by regulator ac- tion. Among the assumptions that were necessary to obtain equation (h.6) was the assumption that the field current remains constant. Kimbark8 and Crary'7 show the development of a formula that alters Eq of equation (h.6) in order to consider these changes in field. The developed formula is: AE E(t+53’-)-0J2L i(t) (1"?) T5: f 0321. i(t)r f + sr f (23) 'r f 2 Eq where the field excitation voltage in armature terms Ef = wfaLsr If L Tf.= fig the field circuit time constant (armature open) f At is the time increment considered Ef is incremented to simulate regulator action and for a given E1, the formula determines the change in Eq or the change in field called the field decrement. For each increment of time, Eq of equation (h.6) is in- cremented by the amount indicated by (h.7) if the change in field effects are to be considered. An example by Kimbark8 uses equation (h.7) to determine the percent change in field for a machine operating with a short circuit at the bus. -37.. The machine parameters for Kimbark's example placed the machine in the same class as the machine of the thesis problem. The results listed for Kimbark's example showed a decrease in field of about 23% in the first half second of sustained short circuit. A third deviation from the direct solution of equation (4.6) occurs if the fault is a less severe fault than a symmetric short circuit at the bus. An example of a less severe fault is illustrated by the problem of a generator connected to a bus with a double circuit and the fault taking the form of a symmetric short circuit on one of the two circuits. In such a situation, equation (h.6) would have different Eq/XdT values for pre- fault steady state, for the interval of the fault, and for the interval after the fault has been removed. Even the equal area criterion.becomes complicated in this third case. That the equal area criterion is not so simple for this type of problem, and that the equal area criterion is not readily adaptable to digital computer solution, is demonstrated.by a qual- itative discussion of an example. Dunne eaeeukt GE“. 7 But tK— Fqut. Figure h.5 The diagram of Fig. (h.5) represents a generator connected to an infinite bus by two parallel lines. A fault occurs on one of the lines. If the faulted line is removed by circuit breaker action, the critical -38- switching angle can be determined using the equal area criterion as fol- lows: Betcvd Afif A FuAL/t g ”rift/ety P au ernnvd’ 03 /—\A,, /7 Davin A. H2 43 Faul 5 5 71' ° c 5 —-> Figure h.6 To apply the equal area criterion here, three different power angle curves are required as shown in Fig. (h.6). The three power angle curves are curves representing the relations in equation (h.3) for the Xd before, during, and after the fault. The horizontal line represents the value of power being supplied to the generator by the prime mover at the initiation of the fault. The cross-hatched areas above the initial power line are taken as positive, and the cross-hatched areas below the initial power line are taken as negative. The angle, 50, is the initial torque angle and 5c is the torque angle at the time the fault is removed. By the equal area criterion, 5c is the critical switching angle when ‘ A1 + A2 + A3 + Ah = 0 The critical 5c is usually found by graphical trial and error means. The switching time is then found by numerical solution of equation (h.6) with Eq and xdT determined for the faulted condition. -39.. 7,8 Pre-calculated swing curves are sometimes used to determine the critical switching time. Pre-calculated curves cannot be applied to the open circuit and short circuit synchronous machines, however. h.h COMPARISON OF RESULTS AND METHODS The program was run for switching times of 0.1, 0.2, 0.25 and 0.3 seconds. Graph h.l of the results indicates that the system was stable for switching at 0.25 seconds and the system was unstable for switching at 0.3 seconds. These results compare favorably with the 0.28 critical switching time obtained by the equal area criterion. In contrast with Kimbark's8 results of a change in field of about 23% in the first 0.5 seconds of sustained short circuit, the if1(t) ob- tained by direct solution here had a maximum variation of about 0.1% of the steady field current If. This low order of magnitude for the change in the field current is as should be expected on the basis of the terms of equations of (h.2). In per unit quantities, the armature and field inductances are: Lf = 5 L5‘37“! and in per unit quantities the voltage disturbances causing the changes in the field and armature currents are of the same order of magnitude. Both are approximately [vd(t) - s L: iq(t)] and did(t) ‘ + 1 __= v(t)-eLi(t) '— dt [ d s q ] L; diflhs) . + x T... = [vd(t) — 9 L8 iq(t)] L; If -40- where K and L; are the same order of magnitude. Kimbark8 represents the so-called change in field flux by a ficticious voltage component in the armature circuit. This voltage component is de- rived from the same relationships that appear in the equations solved in this thesis. In the derivation the voltage component is related to the field current. From the thesis solution, for the time interval of 0.5 seconds, the speed had increased by h.5$. The commonly neglected change in speed appears to be a more significant variation than the field current decrement. Very little more can be said in comparison of the results, since stan- dard techniques usually give only 6 as a function of time, but the id(t) and the iq(t) of the computer solution shown on Graphs h.2 and h.3 exhibit characteristics distinctly different from.the form commonly assumed for the id(t) and iq(t) variables. The approximate mathematical expression for these current variables obtained from the graph are: 1q(t) ids) 0.82 sin (mt + V1) 1.25 + 0.82 sin (wt + Y2) A closer examination of the data showed that the «1 terms were not con- stant at 377 radians/second, but were more nearly 5, the instantaneous angular velocity. It is commonly assumed that id(t) and iq(t) have a slow variation with a period similar to the period of the mechanical oscilla- tions. The period displayed here is around 0.0166 seconds, compared with a mechanical oscillation of around 1.0 second. It is common to assume that for the short circuit condition, iq(t) goes to zero immediately. The amp- litude of the oscillations of iq(t) displayed in Graph h.2 is around twice -hl- the initial value of iq(t). The currents were obtained for the condi- tions that the armature resistance was neglected and currents appear to have the form of sustained oscillations, but the armature circuit time constant is around one second so the decay would be slight in the time shown if the resistance had been considered. A review of the mathematical relations for the equations solved sub- stantiates the results obtained. For the equations (2.6)a and (3.1), with the assumption of constant speed and constant field current, the usual form of Park's equations are obtained: van) -- (rags; £3) 1am +bL;1q(t> (h.8) vq(t) = (112%; £13) iq(t) - 6L: id(t) - bJ—str If Neglect R; and for short circuit conditions, vd(t) = vq(t) = 0, so the equations (h.8) become: 7%;- id(t) = - éiq(t) (as) d , , ET; iq(t) = + e id(t) + 9421.8r If For the initial conditions and coefficients of the present problem, a Laplace transform solution of the two equations yields: iq(t) 0.82 sin (at + 4’1) (n.10) id(t) 1.25 + 0.82 sin (bt + V2) With the resistance considered 0.82 e'o'8t 1q(t) ids) sin (ét + \V1) (”'11) -0.8t 1.25 + 0.82 e sin (8t + ‘V2) It is interesting to note that for the constant speed solution, the 0.82 coefficients of the sine terms and the 1.25 constant are independent -ug- of Speed. .A conclusion that for small changes in speed the id(t) and iq(t) variables would have the same amplitudes and same forms is substan- tiated by the results on Graph h.2. .Admittedly, the average value of iq(t) is zero, and for a high inertia- torque ratio of a power generator, the variation of the iq(t) would add only a slight ripple to the swing curve 6 as a function of time. This high inertia-torque ratio does not always exist, however, when the standard technique of considering a system as an equivalent machine is used. The high-inertia to torque ratio would not exist in small synchronous machines in control systems, so in such a system the standard assumptions would have questionable accuracy. The form.of id(t) and iq(t) is of particular importance when solving the differential equations (2.6)a and (3.1) by numerical methods. If the id(t) and iq(t) variables had a slow variation with a period of the order of one or two seconds, similar to the mechanical period, the time incre- ment for the numerical solution could be taken as a value such as 0.05 seconds or some other value small compared with the period of variation of the variables. The correct form of the currents, id(t) and iq(t), shows a period of around 0.0166 seconds and.a time increment of 0.05 seconds, for the numerical solution could not be expected to give significant re- .sults. Finally, in addition to comparing results of the standard and thesis techniques, consider the differences in principles and methods for the first class of problem. In contrast with the standard technique, the thesis technique: 1. Solves differential equations as differential equations. -h3- 2. Does not require the development or use of the transient reactance concept. 3. By direct solution of equations (2.6)a and (3.1, gives values of armature current, field current, speed, and torque angle as a ftnction of time. The standard method gives the torque angle as a function of time. A. Does not drop field variation terms as in the development of the equation (h.6), so field variation effects can be con- sidered without developing additional formulas. 5. Can solve a large class of problems without changing proce— dure. Changes in specified variables, specified parameters, and switching times are controlled by the master routine of the computer program. V. TWO MACHINE SYSTEMS The question arises as to the desirability of analyzing the single machine system and the two machine system as separate problems rather than as problems included in an n-machine formulation. A partial Justification of this separate consideration is given in the following paragraphs. There are many power system stability problems that fall into the one machine or into the two machine class of problem and give satisfactory correlation between calculated results and observation on the actual sys- tem. The study of the effect of machine parameters on stability, circuit breaker action for connecting or disconnecting individual machines, and many others, are single machine problems. After a preliminary investiga- tion, it may be found that one machine and a system can be adequately rep- resented by a two machine analysis as far as the individual machine is -42.- concerned or two interconnected systems may be adequately represented as two equivalent machines in a two machine analysis. The case actually in- volving only two synchronous machines would arise only in isolated cases or in low power level control systems. The latter is not a power system stability problem. The greater simplicity of the one machine or the two machine analysis methods makes these methods more desirable than a generalized form when these simpler methods give acceptable correlation with observed results from the physical system. The simpler forms occur because of the element- ary patterns for the transmission line systems interconnecting the ma— chines. The multi—machine fonmulation would need to include a network analysis program with the machine analysis program. Whether standard or thesis methods of analysis are used, the two ma- chine problem represents one step in the direction of a multi—machine analysis. Even if methods do not apply directly, information obtained by analyzing a two machine system should be helpful in formulating a multi- machine analysis. Some problems which are conventionally considered as two machine problems are now listed. 1. Two machine representation is used in the study of an individual machine delivering power to a system over a long transmission line. As far as the individual machine variables are concerned, the distant system is adequately represented as an equivalent machine. Studies are commonly made of the effect on the individual machine of a fault and fault clearing time for a fault on the transmission line con- necting the machine to the system. Studies are made of the effect l ‘1}? "W w; -h5- on the individual machine of an increase or decrease of the net gen- eration of the system due to added loads, lightened loads, or faults on the system. 2. Two machine representation is used in the study of two systems which are interconnected by a transmission line. After some prelimi- nary investigation, it may be found that the tie line stability is adequately represented when the two systems are considered as two equivalent machines. Studies are commonly made of the effect on the tie line stability of a fault and fault clearing time for a fault on the interconnecting transmission line. Studies are made of the effect on tie line stability of a net increase or decrease in generation of either system due to changes in load or local faults in either system. 3. The case of two interconnected synchronous machines would be rare in power system study; however, analysis techniques in terms of two machines are commonly presented and serve for the systems repre- sented by two machines. The prOblem worked in this chapter is that of two machines, and under appropriate assumptions does represent two systems. The stability prOb- lem considered is the prOblem caused.by an increment of change in power of one of the machines. The transmission line fault is discussed.but not illustrated by a numerical example. The standard techniques require a different approach to the two types of problems, while the thesis tech- nique is basically the same for the two types of prOblems. 5.1 THE DIFFERENTIAL.EQUATIONS OF'THE TWO.MACHINB SYSTEM In this section, synchronous machine equations are combined, trans- formations are chosen and applied, and the resulting system.of equations for a two machine system are listed. -us- Consider two synchronous machines interconnected by a transmission line system, Equations of the form of (2.1)a and (2.1)b represent the system as follows: (5-1)a V81“) = $831+ a(117-08331 aqtérlwl) ‘Qsl(t) vf1(t) €%”Xr81(¢l) Rfl+ it Lfl if1(t) ’ 4—” of (st) (5-1)b T (t) = - ; [gf (t) i (t) 9 881 srl l 1 2 81 fl ] E73; 4£}81(¢1) Lfl e0 (t) d o 11:“) + (31”13?) “‘1 (5.2). Was) = @882. dc; “a €33.20wa J32“) mm iii-03.32%) Rf2+ a4; 1,1,2 1mm OZ x81‘2(¢2) (5.2)b T2(t) = -% [J;2(t) if2(t)] 5%; 020332 r32 Jags) d o + (B + J -) ¢ 1f1(t) 2 2 dt 2 where the coefficients and matrices are defined as in Chapter II. The (2) subscripts represent variables and coefficients for machine number (2), and, as in the first problem, let the variables and parameters with the (l) subscript represent machine (1) plus transmission system co- efficients’4 and variables. For the specified connection pattern, by use of circuit equations and segregate equations, (5.3) and (5.h) and the application of the symmetrical component transformation of variables, (5.3) (5.4) 32m = 4.200) (2.3), to the stator voltage and the stator current variables, there re- sults (5-5)a where -h7- —va1(ti- —va2(tij ve1(t) vc2(t) —1a1(t3— ‘_iae(t57 :931(t) = :b1(t) = 1b2(t) I. 0 -W Rs+Ls g; 0 + + d 0 0 RS+L8 a"; O = 0 0 vfl(t) o fiLarle‘Jel d - Lyf2(tZJ __ o - EELsrze J92 T1(t) = (Bl + J1 €§>a1 plL srli o o 0 .EL 391 - 1L. 392 ° dthrle dthrae + '* d 11. 'J ._Jl ~J Rs L sdt dtLBrle 61 dtfinfif d e3 EL dtLarleJ 91 Rf1+Lf1 dt 0 d a _.o __ 11(t) + 11(t) 11(t) 1f1(t) 'J 1‘9 m{11(t) e 61} T2“) 3 (Ba * J2 3%)52 " 1’2 Lara 1:: Jm {-1 1“) e .162} R: = R21 + 322 R: = Rgl + R22 L: = L31 + L22 L: = Lgl + L32 9: = wt + 61 + ‘/2 62 = at + 82 + 1/2 and 99:“ {1:(t) {3913 and Jm {41%) {392} means the imaginary part of the bracketed quantities. For the backward sequence transformation of the type represented by equations (2.h), a choice must be made. Some of the forms in the litera- ture that reduce the system to an equivalent machine would seem to sug- gest a transformation in terms of 9, where 6=mt+51+52+1/2 This 6 would be necessary for convenient reduction of the mechanical equa- tions to a single equivalent in terms of 812 = 81 + 52 The overall system of equations was simpler and.the identity of the individual machines was retained, however, when either the 61 or 62 trans- formation was used. Thus, transformations of the form of equations (2.h) were applied to the stator variables for the theta equal to 91' Further, the transformations of the form of equations (2.5) were applied to separ- ate the reals of the backwards sequence variables from the imaginaries of the backwards sequence variables, and the equations (5.6)a and (5.6)b were obtained. o R:+L: §% 0 o o R;+L; §% (5.6)a o = o - élL: vfl 0 (21'er 2%; v o --J§L [cos (5 —s ) ii-- (6 .6 ) sin (a -s ) L f2_ __ 5r2 1 2 dt 1 2 1 2 o o {91 L; JELsrl % R;+L; gt "JEL8r1 b1 0 Rf1+ Lfl dit- . d. ' ' JELsre [5111 (81-62) E + (61-62) cos (51-62)] 0 o ‘ :M r d ' - 2Lsr2 [cos (51-52) at + 62 sin (51-52)] id(t) faLsra [- sin (81—52) % + .92 cos (61-82)] iq(t) o if1(t) Rf2 + L1’2 gt __ _ff2(t) a - f2 . Tl(t) = (51+J1 Et)¢1 + /2 Lsrl p1 1f1(t) iq(t) (5.6)b d . 4; T2“) = (132w2 as, + /2 1:er p irate) [iq(t) cos (61-62) - / ‘ id(t) sin (51-52)] Express the swing equations of (5.6)b in terms of power . 2 (5 7) P = Gi‘JELsrl if1(t) iq(t) + M1 jig 91 = - 52 JELsra 11,2“) [rqm cos (51-52) - 1am sin (51-52)] +M 6 2 dta 2 where d2 on 91=wt+61+x/2 $591 '51 2 — L = o. 92 — wt + 62 + 1/2 dt2 92 52 Equations (5.6)a and (5.7) are the differential equations of a two machine system. The thesis technique is the solution of these, (5.6)a and (5.7), directly by numerical methods on the digital computer. The standard technique makes some simplifying assumptions and uses a modified steady state method to obtain a solution to these equations. 5.2 STEADY STATE INITIAL CONDITIONS FOR THE TWO MACHINE PROBLEM The initial conditions are usually specified as constant speed, con- stant direct current field, constant amplitude sinusoidal armature voltages, a constant average power, and constant power factor. If these initial con- ditions are imposed on the sets of equations (5.6)a and (5.7), the equa- tions (5.6)a are independent of the set (5.7). From the results of the transformations on the steady state sinusoidal armature currents, id(t) = - |Iml sin (51 + B) lq(t) = lIml cos (51 + B) where B is the power factor angle. The results of the steady state solution to (5.6)a are commonly pic- tured on a vector diagram as in Fig. (5.1). l m _ ,l l m Figure 5.1 From the solution of (5.6)a the initial values are: id = - 0.h8h0 51 = 23.6° iq = 0.3350 52 = -280 if1 = 1.5 if2 = 0.852 From.the constant speed assumption: 51 = 377 radians/second éé = 377 radians/second 51 + x/2 = 91 = 1.9823 radians 52 + r/a = 92 = 1.0821 radians O.h116 radians -O.h886 radians The initial values Obtained from (5.6)a satisfy the equations (5.7). 5.3 NUMERICAL STUDY OF THE STABILITY PROBLEM OF A.TWO MACHINE SYSTEM In this section the constant field current assumption is Justified. The equations are put into standard form for numerical solution. The -52- computer program is discussed, and the use of the program is illustrated in solving the second problem. Upon examination of the results of the first problem, it is found that the percent variation in the field current was less than 0.1 percent of the steady direct current. The system inertia was close to a.minimum, and the fault one of the most severe likely to be encountered. .As these results suggest, examination of terms in the differential equations shows certain terms to be negligible. Some of these terms will be discussed and eliminated as the two machine differential equations are put into standard form for solution. l 0 +4.; Lfl - 2 ETCOS(51-52) 0 1 0 -42 gains -5 ) rt 1 2 $21.31.1 0 l 0 Lfl JEL L sin (5 -5 ) 81‘ 81‘ (5.13) 2 eos(o -52) f2— 2 1 2 0 l Lfe 1 Lf2 0 0 0 o 0 o 0 0 o 0 0 0 __ 0 0 o 0 —f "" 0 0 o 0 dt 1 C1(t) 0 0 0 0 ad; 1q(t) 0 0 0 0 34,—; 1:1(t) 0 0 0 o 5% if2(t) = 1 0 0 0 '91 0 1 0 0 "l 0 0 l 0 '92 0 0 0 14 _ 92 J -53.. + u + O l [-1281 d - elLsiq+ (21451-292 sin (51—52) ire] L; T + S a + o ' 1 [-R 1(1 + elLsid+ ~f2Lsrlel ifl - ~[2L3r2 92 cos (51-52) ire] I: l [Vfl ‘ Rfl in] Q; - J2Lsr2(él-52) sin ($1.52) 1d - J2Lsr2(él-52) cos (ol-oanqwfz-Rfeifa] it .91 [P1 - bl J21.“1 ifl iq] “11 2’2 ° 1 {P2 + 92 J2L5r2 if2 [iq cos (51-82) - id sin (81-52)]; E L— The 8 x 8 coefficient matrix of the first derivatives is of the form 6 o 11 0 B 22 where the 5 is a unit matrix. 22 The inverse is obtained when 61-: is obtained. The submatrix 611 is a 1|» x h and the inverse will not be presented here but significant terms will be discussed. The determinant of 811 is 2 2 2 2 Lara 2 Lsrl + (2Lsr1 L81?) sin2 (8 -5 ) " + ' + +2 1 2 L12 Ls Lfl Ls Ls Lfl L1’2 but _ l J—2Lsr1 = JELsrg - 3T7 p.u. L+ = 3— . . s 377 p u and the order of mgnitude of the field inductances would be about 6 p.u. for each Lfl and 1’12“ Thus, the determinant of 611 is approximately unity. ~5h— Similarly, all the diagonal co—factor terms of the inverse are approx- imately unity. For illustration purposes, consider the 33 co-factor a33. 2 12 2L 2(-——~) 833: __s_rf =1__317.2_ =1-6——1§—7? =1-.00088=1 _— x LfaLs 6x377 The off diagonal co—factor terms a12,a13, and 81% which reflect the right hand term into thed -— ti d(t) equation are: 2 — 2Lsrai(oo) (so) 312 - Lf L: 3n 1" 2 cos 1" 2 2 s 2 . libs; - 1m sin-egg) 13 L+ L L+2 l 2 f2 5 EL a = sr2 cos (5 -6 ) 1h + 1 2 B The first four lines of the right hand side of (5.13) are of the form [- ‘1 A D..- _ L1?) where B, C, D are of the sane order of magnitude as A or smaller than A. Thus, (t) =A/L: + 13‘;/L + c/L + D/L dt 1d 312 8‘13 f1 311. 1’2 where the last three terms are approximately (.OOOhh)(A/L;). A similar analysis applied to the a; iq(t) equation leads to similar results, so the equation (5.13) in standard form for solution with field variables as constants is as follows: '55- ~— — ”' '1 d + 0 + o , l Etid [’Rsid‘elLs i{6’2 ~[xi-"1:2 “1“ (51'52) 1:2] L;- d + . + ° - 1 atiq [-Rfiiqi-elL8 id+91 J2L3r1 if1'92 J— 21"er cos (51-52)if2] L: 6‘1 91 (5.1%) = co 0 . 1 91 [P1 " 91 JZLBrl if). lq] If; 92 e2 .9. P+é~f2L i icos(5-5)-isin(8-8)]-1=— L_2_J 228r2f2q 12 d 12142 The equations (5.1h) are the differential equations that describe the two machine system. .Equations (5.1h) are programmed as the closed fr sub- routine for the F-6 differential equation solving routine. The F-6 and the ff routine for (5.1%) under the control of a master routine can be used to solve all of the two machine stability problems discussed in this chapter. The master routine as in the first prOblem controls the amount of information extracted from.the computer, simulates faults and fault re- moval by changing parameters in the fr routine, and simulates changes in load by changing constants. For the particular prOblem, the master routine was written to read out t, 51, 52, id and iq. The change in load was affected on read in by determining initial conditions for P1 = - Pé = 0.5 p.u. power and reading in that P2 = - 0.6 p.u. at the beginning of the program. The frequency of the read out and the duration of the program are controlled by the master routine. With the f} routine written for a two machine system, and a master routine written for the type of stability prOblem.under consideration, the second problem is solved.when.the proper initial conditions and system parameters are supplied.to the computer with the program. -56. For the second problem, the scaled initial conditions and parameters are supplied to the computer as follows: The numbers in memory location 3 through 8 as required by F-6: In 3 In h In 5 In 6 In 7 In 8 In 10 The initial The integer 20 indicates the location of the first of the sequence of initial conditions. The integer 27 indicates the sum of entry 3 and n, the number of equations in the system. The integer 13 indicates that the scale factor on the fr is 2'13. The integer 10 indicates that the time increment is 2-10 see. The integer 36; The g? will have their error held to about 1 in the 361"-12 binary bit. The integer 160 gives the memory location for the entry into the fr routine. The integer 100 indicates the memory location of the first of a sequence of scaled coefficients and constants. values of all the variables are scaled‘by 2"13 and read in the following location. Location 20 21 22 23 2h 25 26 Variable Scaled value t O 52 0.01.602 92 0.0001322 b1 0.0h602 91 0.0002h20 iq 0.0000h088 i - 0.0000590 -57... The coefficients and constants for the fr routine: Location 100 101 102 103 10h 105 106 107 108 109 110 111 112 113 Quantity Pl/Ml Fe + AP2 Scale Scaled Value 0.001917 0.006903 0.0000502h 0.00005965 0.001h6h8 0.0008320 0.03183 0.01061 0.01 0.36816h 0.73632 0.0001917 0.63662 0-005305 For the master routine, supply the following counters to control read out. If the integer "k" is placed in the 301:11 word of the master routine, read out will occur after each (k+1)st calculation. If the integer "p" is placed in the 36th word of the master routine, read out will occur."p" times before the programistops. The pro-set parameters in memory locations 3 through 8 must be read in before F-6 is read in. end of the program tape. All the other information can'be placed at the At the end of the read in, transfer control to the left side of the first word of the master routine and the computer I .-'§ 0 I 3 s A ' -' -Mfiat‘sfl‘, m.- {u ‘ ' -58.. gives the results. A common stability problem is the problem of determining the maximum AP2 that can be permitted and still retain stability. The maximum AP2 can be determined by running the program a number of times with different values of P2 + AP2 in memory position 101. The assigned values of AP2 can be mde to converge on the critical value. The problem was run here for AP2 as 0.1 p.u. and 0.3 p.u. or P2 + A132 as 0.6 and 0.8 p.u. power. The results are presented on Graphs (5.1), (5.2), (5.3) and (5.1;). 5.1; STANDARD SOLUTION Fm THE TWO MACHINE PROBLEM In the literature there are various approaches to the two machine problem. The approach varies somewhat with the particular problem. If it is desirable to consider system losses, the formulation is different from the formulation for the lossless case. For determining critical switching times, equal area criterion or pre-calculated swing curves may be used. For stability information with added load increments, a form of the equal area criterion may be used. All the standard techniques are based on the concept of transient reactances and voltages behind transient reactance , however. To obtain the voltages behind the transient reactances for the two machines, the steady state initial condition solution of section (5.2) is a necessary first step. For the second step, each machine is considered separately as in the case of individual nachines of Chapter IV. From the initial condition, 1 d’ i , vd, vq, and machine transient reactances, a Q voltage behind transient reactance is found for each nachine with equation (hJ-t). These voltages are designated here as El and E2 for mchine (1) and machine (2) respectively. Next, if the mchine and transmission line all .fln l‘ i -59- apnoea“. S N .’ $9360.“. Q.‘ w». J, -60- \.\ QV {attest kx N m6 ”6 m6 9% MG V6 “Q N6 \6 0.0 “In! 01 fun/n) ‘ p.utouuh. S N . . m6 Wm mm “.6 N6 \6 0.6 -62- 770’ u/ ,wzJJng $6 36 36 mic N50 $2303... S V. \\.Q 9.0 m8 86 NQG 86 fig 3Q MQQ N66 \8 '7770’ u/ 37' -6h- losses are neglected, the system can be reduced to an equivalent single machine form. 2 d 5 IE1||E 2| l2 __ (5.8) Mo 1? - Pl“ x12 Sin 512 where L2 M =M1 + M2 6 is the angle between E and E initially, and the angle 12 1 2 between machine rotors during the transient. X12 is the sum of the transient reactances between points, where E1 and 32 are theoretically measured. Equation (5.8) has the same form as the swing equation for a single machine on an infinite bus. Under the assumptions that were necessary to obtain (5.8) , the two machine problem can be solved by the standard. method discussed in Chapter IV for the type of fault that would alter the tran- sient reactance. If 512 increases without limit, then the machines fall out of step. Equal area criterion can be used. If losses are to be considered, the equivalent machine form of equa- tion (5.8) cannot be obtained. A somewhat more complicated form results. 2 d 5 P - T P - T 12 _ 1 e1 __ 2 e2 (5'9) 2 " [ ] dt “1 M: where 2 IEll IEillfial Te1 | 3 11' sin “11* W sin (812 4:12) (5.10) E212 IEl_____l____IE2| Te2 = '3 22‘ sin 0122 ' ,3 12' sin (812 + (112) and where -65- ~911 and ;?22 are driving point impedances (transient) $212 a transfer impedance (transient) all’ 0:22, 012 are the complements of the respective impedance angles. For problems in determining critical switching times, the equal area criterion can be applied using (5.9),‘but the power angle curve is not a sinusoid and the graphical plot and analysis is tedious and subject to inaccuracies of graphical methods. If losses are neglected, the equiva- lent machine form of equation (5.8) is Obtained and the critical switching time prOblem is the same as for a machine on an infinite bus if the equal area criterion is used. For the particular problem of this chapter, the equations that per- mit consideration of changes in generation or load for machine (2) with losses considered, are: 2 d 812 _ Pl - Tel P2 + AP2 - Te2 (5.11) 2 .. ——-—- + M dt M1 2 If losses are neglected P1 = ' P2 Te1 = -Te2 then a? 812 AP2 Ml + M2 (5'12) :5- = “a; + ‘12??? (P1 - Tel) = G(8’12) With losses, the right hand side of (5.12) is not a simple sinusoid as in the case or the individual machine on an infinite bus, but a plot of the right hand side of (5.11) as a function of 812 gives a figure as in Fig. (5.2), where the so is the initial value of 512. -66- A I 63(5/m> ‘9 . As I// ./ / / 50 ////////’ 5);“) Figure 5.2 By the equal area criterion the system is stable if the increment of load is such that A1 < A2. For the machine of problem two neglecting losses: H1 = 6 H2 = 2 AP2 = -0.1 .r 0(512) = [+ Qé—l + % (0.15 - 1% sin 512] = [- 0.05 + .333 - A26 sin 512] = .383 - .h26 sin 612 But at 0(812) = 0 at 612 = 6h° 6&0 A1 = 51 6 [.383 — .1126 sin 512] d 512 = .0072 = .026 180-6!+° A2 = [61.0 [.283 - A26 sin 512] d 012 Thus, for the increment of load AP2 - -O.1 and the lossless assumptions, I A 2’1‘ E . '_ fl: " - . ~. » ’ to ”2'. -. ‘m‘e’c'... ' , ‘_ 5‘“. fl ‘- . ghzg“.. ...' 4‘. w—w—w ,-‘ v '1 .-' -67- the equal area criterion indicate stability. For the increment of load AP = - 0.30 2 1f 0(512) = 1:9:§.+ {% (0.5 - 1‘5 2 '852 sin 51;] at 0(512) = [.u83 - .u26 sin 512] but the curve does not cross the 1f G(512) = 0 axis, so A1 >’A2 and the system would.be unstable for AP2 = - 0.30 p.u. power. A common problem in stability studies is the problem of determining the maximum allowable AP2 for a given initial steady state power fer sta- bility to be retained. A continuation of the above calculations would permit convergence on such a critical AP2. If losses are considered, the solution for the point at which the curve crosses the axis and the evaluation of areas would have to be per- formed graphically. 5.5 COMPARISON OF RESULTS AND METHODS Graph (5.1) of the results, a plot of 8 - s 1 2 illustrates that the system is stable for the case where AP: is 0.1 p.u. and that the system is unstable for the case where AP2 is 0.3 p.u. These as a function of time, results correspond with the results of the standard solution. Since the standard solution to a problem of this type simply indicates if a system is stable or unstable for a particular AP2 and does not give information on 81, 52, 51 - 52, id, and iq, no further comparison can be made. The additional information available from the thesis solution war- rants some discussion, however. The fact that 51, 8 and.6 - 6 are given as a function of time in 2 l 2 the thesis solution, contributes information on how the individual machine -68- or system contributes to the stability characteristic of the combination. In the standard approach this information is not available if the various criteria are used. When the individual 8's are desired, the standard techniques require the multi-machine approach and a numerical solution to two simultaneous swing equations of the form of equations (5.10). Also evident from the Graph (5.2) is the amount'both machines fall behind the synchronous position even when the system is stable. As in the first problem, the current variables id(t) and iq(t) have a form much different than the form usually assumed for these variables. Graph (5.3) is a plot of the envelope of the id(t) and iq(t) variation for both AP2 values considered. Graph (5.h) is a plot which pictures the variation of iq(t) for the case of AP2 = 0.1. The oscillatory component of the current variables is not so pronounced as it was for the circuit fault, but the variation is large enough to be an important factor in choosing the program increment. In any stability problem, an important consideration is the current variation in relation to the circuit breaker action or relay action. Stand- ard solution for the two machine problem requires a separate solution for the currents when relay action is being considered. The solution presented here gives the currents in the direct and quadrature axis component forms. If desired, an inverse transformation to convert these to terminal vari- ables could be included in the master routine. In noting the differences in principles in standard solution and the solution presented here for the second problem, these differences are practically a repeat of those listed for the first problem. In contrast with the standard techniques, the thesis techniques: l. Solves differential equations as differential equations. 2. Does not require the development or use of the transient reactance concept. 3. By direct solution of equation (5.1h), gives values of armature current, 8 and 82, as functions of time. 1; h. Can solve a large class of problems without changing pro- cedure. The programmed equations (5.1h) represent the system. Different stability problems are simulated by using the proper master routine. 5. Considers losses or neglects losses without a change in basic form. VI. SUMMARY The common simplifying assumptions were considered in two groups in this thesis. The first group of simplifying assumptions were the assump- tions that fundamental harmonic representation is accepted and that the resistance and inductance coefficients can be represented by constants. For the first problem with the first group of assumptions, the synchronous machine on an infinite bus is represented by the non-linear differential equations (2.6)a and (3.1). For the second problem with the first group of assumptions, the two machine system is represented by the non-linear differential equations (5.6)a and (5.7). The technique of this thesis consists of solving these systems of non-linear differential equations with the use of the digital computer and numerical analysis methods. A second group of simplifying assumptions were necessary in order to obtain formulas used in standard stability techniques. These assumptions were constant speed, constant direct current field, and steady state armature .r- ’31”: " -70- voltages and currents. The thesis solution is more general, in that the second group of assumptions are not necessary, and further, the identity of all the variables is maintained for more detailed study. For the electromechanical transient of a stability study, the numeri- cal solution of the equations (2.6)a, (3.1), (5.6)a and (5.7), gives the stability information directly. The standard solution of the electrome- chanical transient problem utilizes the assumptions that the electrical variables can be described by a steady state algebraic system of equations with the coefficients altered to values which apply under transient condi- tions. These transient parameters and fictitious armature voltage compo- nents are used in the numerical solution of the mechanical equation or the swing equation. The results of the computer solutions illustrated that of all the second group of assumptions, the assumption of the steady state character of the electrical variables was most radically out-of-line. Standard solutions of stability problems often use criteria based on the second group of assumptions and the further assumption that the system losses can be neglected. In contrast with the solutions to similar problems in the literature, the field current variation was shown to be negligible in the results of the first problem. It can be concluded that the assumption of constant field current is a reasonable assumption; however, rather than eliminate the field current term as is done in standard solutions, it is convenient to leave the field current terms in the equations so that these terms can be incremented to simulate regulator action. The change in speed was found to be a more pronounced variation than the variation in field current. The constant speed assumption did not -71- appear to be unreasonable even for the low inertia and severe fault of the first problem, however. The fact that the thesis solution gives the additional information about the field current, speed, and armature current components, while the standard solution does not give any of these directly, is certainly significant. The armature currents must be determined for relay and cir- cuit breaker study. The variation of iq(t) and id(t) illustrated in the results of both problems was probably the most pronounced deviation from accepted theory. Clearly, both id(t) and iq(t) have forms or components containing the fonm of II] an(ét+ y) In the standard derivation and use of the id(t) and iq(t), it is common to start with Park's equations, _ + + £L_ - + vd(t) — (RS+L3 at) id(t) + 9 L8 iq(t) (6'1) + +d ' + 'J- vq(t) = (R8+Ls a?) iq(t) - 9 L8 id(t) - e 2LBr If Certainly it is not reasonable to assume the terms + d (6 2) Ls dt 1d(t) + d Ls at iq(t) are zero or negligible for the current variations displayed in the results presented here. It can be seen, however, that neglecting the relations (6.2) in (6.1) and redefining ~ + e L i (t) (6.3) . 2 ‘1 e 1.8 id(t) -72- in terms of transient reactances could lead to results similar to those obtained by direct solution of (6.1). The form of the current variables played an important part in.the determination of the proper time increment for the numerical solution. In particular in the first prOblem, the time increment had to be small compared with the period of the current variation, approximately 0.0166 seconds. The standard techniques for both the one machine and the two machine problems require that different methods or criteria be used for different faults or different stability problems. If losses are to be considered, still a different form is needed. In sharp contrast, the thesis procedure is very flexible and is adapted readily to many types of stability'prob-' lem. The ff routine is programmed as required by the particular system under consideration. There are trivial differences in programming for the lossy and the lossless cases. 'With the system adequately described by programming f}, all the various types of faults and stability prOblems can be solved by using a pr0per master routine. The methods of solution used here and.the advantages of these methods illustrated.here extend to areas other than the areas of power system.sta- bility. The synchronous motor connected to an infinite bus has the same differential equations as the generator on the infinite bus. A.master routine that controlled the Pin term could use the one machine program to study the synchronous motor for various periodic or impact loads. The control system problem of two synchros has the same differential equa- tions as the two machine system, The terms sometimes neglected in power systems may not be negligible with control system.apparatus, however. -73.. The equations and programs presented here are for round rotor ma- chines and for symmetric faults. The programs also cover only one and two machine systems. Further investigation aimed toward extending these methods into the areas of salient pole machines, non- symmetric faults, or multi-machine systems, appears to be indicated by the results pre- sented here. .A conclusion based on the results of the investigation reported here is that, for the class of problems considered in this thesis, methods de- veloped before the advent of the digital computer should not be used as the basis for digital computer investigation. The more fundamental syn- chronous machine differential equations can now be solved directly. APPENDIX A An Existence Theorem E. L. Ince5 gives an existence theorem.which applies to the problem considered here. Appendix.A is a restatement of that theorem. The system of equations to which the existence theorem is to apply is a system of n ordinary differential equations in the n+1 variables. x, yl’ ya, 0000000, yn0 The desired solutions are the n equalities 3’1 = 81(X); 3’2 = 820:). H": yn = sub!)- These equalities are to satisfy the conditions 0 _ o o _ o o = o yl- 3].(x )’ y2— 82(x )’ 000000000, yn %(X) where x0 0 o o ’ yl’ ya, 000000000, yn represent the initial values of the variables. Further, for the theorem to apply, the differential equations must be expressible in the normal form: 1.. ’53:" £10.. 1'1. ya. ......, yn) “'2 “at; " f2(x: 3'1: 3'2) 00-0-0: yn) (11.1) na— 75?" fnh‘: 3'1) 3'2: ......, yn) Equations involving higher order derivatives than first order are included by’a change of variables. Specifically, consider, for example, the equation 2 d _ E (21.2) 3% - fix, an iv) -A.l_ To obtain the required normal form, let 2-91 dx and write (3.2) as a pair of simultaneous equations. (A.3) Theorem = f(x: 3') z) EH? E13 a. .1: For a system of equations of the type of (3.1), let (x0, yi, yg, ......, y:) be a set of real numbers assigned to the real variables (x, yl, ya, ......, yn). Let D be a domain defined by the inequalities: o o o IX‘X Isa, ’yl-yllsbl’ COD-0000’ lyn-ynlsbn- Let M be the greatest of the upper bounds of fl’ f2, ......, fn for arguments restricted to D, and let h be the least of a, bl/M, b2/M, ......, bn/M. If f1 (n+1) arguments when the arguments are restricted to D and if , f2, ......, fn are single valued and continuous in the the Lipschitz conditions lfr(x, yl’ y2, 00000, yn) - f(x’ Y1, Y2, 0.0000, Yn) I < 1(1- lYl-yll + K2 IY2-y2l + sconce-o + Kn Yn-ynl apply for each r = l, 2, 3, ......., n when (x, Y1, Y2, ....., In) and (x, yl, 32: ......, yn) are in D, then there exists a set of unique continuous solutions, vi = 21(X), V2 = 32(X). ......., yn = sn(X)o These solutions satisfy the differential equations for all x such that Ix — xol < h, and reduce to: -A.2- II N O O 0 3'1, ya, 000000, yn for X The types of terms occurring in the synchronous machine equations (5.7) and (6.11;) are included in the following: f = A1 sin kx + A2ylyh + A3y2 1 f2 = Bl cos kx + Bayuy3 + B3yhy2 + Bhyl f3 = Cl sin kx + Cayhyi f5 = El + Eey'hy3y1 Clearly, the elementary functional forms in (3.h), are single valued and continuous functions for finite (x0, yi, ....., y?) and for finite (a, b1, b2, ...., b5) defining the domain D of the theorem. Further, for these types of equations the Lipschitz conditions apply for the arguments in D. The domain, D, of the existence theorem is, in general, very small relative to the range of the variables to be considered in any given prOb- lem. Different initial values of the variables and corresponding domains covering the range of the variables extend the existence theorem to a larger range of the variables. Even with the assurance that a unique continuous solution exists for a desired range of variables, numerical approximation solutions are subject to limitations and errors of the particular numerical method used. APPENDDi B The Operation of the Digital Computer Program The digital computer library routine, F-6, uses the fourth order Runge-Kutta formulas. The Runge-Kutta formulas, F-6, the differential equation routine, and the master are described in Appendix B. For the system of differential equations in the normal form of equa- tions (A.l), the r32 equation is the first derivative of the rth variable: dt y = fr (t) 3'1: _-'—-) y 1- , yk) r For given values for all the variables at the beginning of the nth incre- ment of time, the Runge-Kutta formulas determine the value of the variable yr at the end of the nth increment or the beginning of the (n+1)at incre- ment. Let the subscripts indicate the variable, the superscript the inter- val, and fr the first derivative of the rth variable, then by Runge-Kutta (n+1) r n (B.l) yr — yr + Ayr where = 1 I _1; II I]; III 1 1v (B.2) Ayr B Ayr+3 Ayr+3 A yr+6 A yr and o n n A yr = fr [15 ) V1) """ 2 Y?) "") YR] At .. At A__’y1 A'y n 43% A yr = fr [tn + 2, y§+ 2 —, ---, yr + 2r, ---, yk + 2 ] At (B-3) u n at A y A y A at n n l n r n 3k Ayr=fr[t+ %'*1y1 2’---’yr+ 2’---’yk+ 2]“? Aivyrf = frPt + At. Y1 + A” yl. ---, y: + AmyJr . --- yk + A” y k] M The fr for this thesis problem are the expressions of equation (h.2). The digital computer library routine, F-6, used here controls the amount of truncation error and the increment by comparing Ayf of expression (B.2) with A"yr of equations (B.3). If these two approximations, Ayf and b.\ -B.2- A"'yr, for the increment of the variable yr differ by more, for any yr, than an amount specified by the user, the program recomputes Ayr with the interval reduced by half. If necessary, the routine continues de- creasing the increment until the difference is within the specified limits. The user of the F-6 routine must supply a closed subroutine which evaluates the fr of formula (A.l). The F-6 routine referes to the fr routine to evaluate A'yr, A"y£, A”’y and Aivyr. Scaling must be used r’ so that all fr remain less than one for the range of the variables con- sidered. In order to use F-6, the digital computer memory locations 3 through 8 must contain the following parameters during read in and operation of F—6. Location 3 must contain the number "a", where "a" is the memory lo- cation of the first of the sequence of locations for the initial values of the variables. The t initial value is in "a" and the initial values of the other variables follow in sequence. Location u must contain (a+n) where n is the number of equations in- cluding dt/dt = 1. Location 5 must contain the number "m". "m" is an integer and the fr routine calculates scaled fr’Er’ such that fr = 2'“1 fr < 1 location 6 must contain the integer to - m. The quantity 2‘30 is the specified scaled increment At. The program decreases this increment if necessary. Location 7 must contain the integer "e" which is a number such that the integral part of 3/h e is equal to or less than (50 - m). The integer -B.3- "e" is used to specify the required accuracy of the calculations. Location 8 must contain the integer "b" which is the memory location of the entry into the fr routine. When the master routine transfers control to F-6, F-6 uses the scaled initial values of the variables and the fr routine to compute new scaled values of the variables. The initial values of the variables are replaced by the new values, then control is transferred back to the master. The fi routine is the program.of the differential equations (A.l). The differential equations represent the system, thus the ff adapts the program to the system. The fr routine must be written as a subroutine, which, when entering, evaluates the scaledfr and then returns control to F-6. The fr routine takes the scaled variables from the computer memory locations containing the initial values required by the F-6 program.and calculates the scaled fr according to formulas (All). These fr are placed in the memory positions a + n + r, where "a" through "a + n - 1" contain the initial values of the variables. The evaluation of the fr for this particular problem requires the determination of vd(t) and vq(t) where vd(t) IVI sin.8 vq(t) - IVI cos 6 In order to evaluate the trigonometric functions, the f} program uses the digital computer library routine for sine and cosine evaluation. The fr program was written so that it could be applied to any machine operating on an infinite bus. The integer "g" must be placed in the digital com- puter memory location 10. "g" is the location of the first of the scaled coefficients and constants of formula (A.l). The other coefficients are scaled and read in in a definite sequence. The scalings and the sequence —B.h- are presented in the solution of the problem. After each increment of time, F—6 returns control to the master rou- tine. Thus, the master can be written so that it removes the scale fac- tor, prints any of the desired variables, and prints out these variables as frequently as would be useful. At any desired time, the time deter- mined by counting the increments, the master routine can simulate faults, circuit breaker action, or changes in load conditions. The master simu- lates these actions by changing coefficients or constants in the fr rou- tine at a Specified time. For example, the short circuit at the bus in the first thesis problem is represented by the bus voltage being set equal to zero. At a time, determined by a counter in the master routine, the bus voltage is returned to the pre-fault value. The action of returning the bus voltage simulates circuit breaker action and the action is per- formed by the master routine. 2. 10. ll. 13. BIBLIOGRAPHY Blandel, A. , "Transactions of the St. Louis Electrical Congress," 190h. Park, R. H. , "Two-Reaction Theory of Synchronous Machines - Part I, Generalized Method of Amlysis," A.I.E.E. Transactions, Vol. 1+8, pp. 715-30; Jlfl-Y: 1929- Koenig, H. E. and Blackwell, W. A., "On the General Properties of Terminal Equations for Polyphase Machines," Final Report Project G-3062, National Science Foundation. Koenig, H. E. and Blackwell, W. A., book unpublished. Ince, E. L., Book - "Ordinary Differential Equations," Dover Publica- tions, Inc. Puckstein, Lloyd, and Conrad, Text - "Alternating Current Lhchinery," Wiley. Crary, S. B. , "Power System Stability, " Vol. II, Text, Wiley. Kimbark, E. W., "Power System Stability," Vol. I and III, Text, Wiley. Standard Handbook for Electrical Engineers, A. E. Knowlton, McGraw- Hill. Kunz, K. 3., "Numerical Analysis," McGrav-Hill. Johnson, .D. L. and Ward, J. B., "The Solution of Power System Stability Problems by Means of Digital Comuters," A.I.E.E. Transactions, Vol. 75, Part III, 1956 - pp. 1321-7. Brown, W. T. and Cloves, W. J ., "Combination Load-Flow and Stability Equivalent for Power System Representation on A-C Network Analyzers, " A.I.E.E. Transactions, Vol. 7h, Part III, pp. 782-6, 1955. Narinder, E. 0., "A Rational Method for the Step-by-Step Calculations in Power System Transient Stability Studies," A.I.E.E. Transactions, Vol. 7h, Part III, 1955, pp. 1087-91. lit. 15. 16. 17. -2- Lane, C. M. , Long, R. W. , and Powers, L. N., "Transient Stability Studies II, Automtic Digital Computation," A.I.E.E. Transactions, Vol. 77, Part III, 1958, pp. 1291-5. Gabbard, J. L. , Jr. and Rowe, J. E., "Digital Computation of Induc- tion-Motor Transient Stability, " A.I.E.E. Transactions, Vol. 76, Part III, 1957, pp. 970-5. Shackshaft, G. and Aldred, A. 8., "Effect of Clearing Time on Synchro- nous Machine Transient Stability, " A.I.E.E. Transactions, Vol. 76, Part III, 1957, pp. 633-7. Shankle, D. E., Murphy, c. M., Long, R. w., and Harder, E. 1., "Tran- sient Stability Studies ," Part I - Synchronous and Induction Machines, A.I.E.E. Transactions, Vol. 73, Part III B, 1951‘, pp. 1563—80. Y 5'98 1%}. t J ' ‘ a... "a 8" ’ i all 1 sum-vs i't‘ iv," 5% E a CC 31v S t Jars-s llflllllfllHllWlH 198 11W 1293 03 3