SYNCHRONOUS GENERATOR MODELS FOR THE SIMULATION OF ELECTROMAGNETIC TRANSIENTS by V l a d i m i r Brandwajn B . S c , Technion, I s r a e l i I n s t i t u t e of Technology, 1971 M.Sc, Technion, I s r a e l i I n s t i t u t e of Technology, 1973 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENT FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n the Department of E l e c t r i c a l Engineering We accept t h i s t h e s i s as conforming to the r e q u i r e d standard THE UNIVERSITY OF BRITISH COLUMBIA A p r i l , 1977 c V l a d i m i r Brandwajn, 1977 In presenting th is thes is in p a r t i a l fu l f i lment of the r e q u i r e m e n t s f o r an advanced degree at the Un ivers i ty of B r i t i s h Columbia, I ag ree that the L ibrary shal l make it f ree ly ava i lab le for r e f e r e n c e and s t u d y . I fur ther agree that permission for extensive copying of t h i s t h e s i s for scho la r ly purposes may be granted by the Head o f my Department o r by h is representat ives . It is understood that c o p y i n g o r p u b l i c a t i o n o f th is thes is f o r f inanc ia l gain sha l l not be allowed without my wri t ten permission. Department of z . The Univers i ty of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date 1 ~r\ 1 ABSTRACT Techniques for modelling of synchronous generators in the simulation of electromagnetic transients are described. F i r s t of a l l , an adequate mathematical model of the generator is established. It uses the conventional set of generator data only, which are readily available, but i t i s flexible enough to accommodate additional data, i f and when such become available. The resulting d i f f e r e n t i a l equations of the generator are then transformed into linear algebraic equations, with a time varying coefficient matrix, by using the numerically stable trape-zoidal rule of integration. These equations can be interfaced with the equations of an electromagnetic transients program in one of two ways: (a) Solve the equations of the generator simultaneously with the equations of a three-phase Thevenin equivalent ci r c u i t of the transmission network seen from the generator terminals. (b) Replace the generator model with a modified Thevenin equiva-lent c i r c u i t and solve the network equations with the gener-red ator treated as known voltage sources e ^ (t-At) behind con-red stant resistances [R , ]. After the network solution at each pn time step, the stator quantities are known and used to solve the equations for the rotor windings. These two methods cover, i n principle, a l l possible interfacing techni-ques. They are not tied to the trapezoidal rule of integration, but can be used with any other implicit integration technique. The results obtained with these two techniques are practically identical. Inter-facing by method (b), however, is more general since i t does not re-quire a Thevenin equivalent c i r c u i t of the network seen from the generator i i terminals. The numerical examples used i n this thesis contain compari-sons with f i e l d test results in order to verify the adequacy of the generator model as well as the correctness of the numerical procedures. A short discussion of nonlinear saturation effects i s also presented. A method of including these effects into the model of the generator i s then proposed. Typical applications of the developed numerical procedures include dynamic overvoltages, torsional vibrations of the turbine-generator shaft system, resynchronization of the generator after pole slipping and detailed assessment of generator damping terms in transient s t a b i l i t y simulations. i i i TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS iv LIST OF ILLUSTRATIONS v i LIST OF TABLES v i i i ACKNOWLEDGEMENT ix 1. INTRODUCTION 1 2. SYNCHRONOUS GENERATOR MODEL 3 2.1 General Remarks about Physical Device Modelling 3 2.2 Model of the Electric Part 3 2.3 Calculation of Parameters for the Model of the Elec t r i c Part 12 2.4 Recent Proposals for Improvements in Parameter Accuracy . 17 2.5 Model of the Mechanical Part of the Generator 19 2.6 Conclusions 23 3. NUMERICAL SOLUTION OF THE GENERATOR EQUATIONS 25 3.1 Choice of Integration Method 25 3.2 Physical Interpretation of the Trapezoidal Rule of Integration for a Lumped Inductance 30 3.3 Choice of Coordinate System 36 3.4 Multiphase Equivalent Networks 38 3.5 Three-Phase Equivalent Circuit of the Generator 41 4. INTERFACING THE GENERATOR MODEL WITH THE TRANSIENTS PROGRAM . 45 4.1 Problem Formulation 45 4.2 Method I - Interface by Means of a Thevenin Equivalent Circuit of the Transmission Network 47 4.3 Limitations of Method I 50 4.4 Method II - Interface with a Modified Thevenin Equivalent Circuit of the Generator 51 4.5 Remarks about Method II 56 4.6 Numerical Examples 60 iv Page 5. INITIAL CONDITIONS, DATA SCALING AND SATURATION 73 5.1 C a l c u l a t i o n of the I n i t i a l Conditions f o r a Synchronous Generator 73 5.2 Consistent Per U n i t (p.u.) System and Conversion to P h y s i c a l U n i t s 77 5.3 S a t u r a t i o n i n the Steady-State Operation of a Synchronous Generator 83 5.4 D e f i n i t i o n s of S a t u r a t i o n i n the Simulation of E l e c t r o -magnetic Transients 86 5.5 Implementation i n the Transients Program 90 6. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH 92 REFERENCES 94 APPENDIX 1 99 APPENDIX 2 103 APPENDIX'3 \ . . 105 APPENDIX 4 107 APPENDIX 5 109 APPENDIX 6 112 APPENDIX 7 116 APPENDIX 8 118 v LIST OF ILLUSTRATIONS Figure Page 1 Schematic r e p r e s e n t a t i o n of a synchronous generator. 5 2 Comparison of the simulated f i e l d current i ^ . 16 3 I d e n t i c a l r e s u l t s f o r current i c i n the f a u l t e d phase 17 w i t h approximate and accurate parameter conversion. 4 T o r s i o n a l model of a turbine-generator u n i t . 19 5 Schematic r e p r e s e n t a t i o n of a three r o t a t i n g masses 25 system. 6 Three phase-to-ground f a u l t at generator t e r m i n a l s . 28 7 Comparison of s i m u l a t i o n r e s u l t s f o r the f i e l d current 29 i n case of a three-phase f a u l t . 8 Comparison of s i m u l a t i o n r e s u l t s f o r the f i e l d current 30 i n case of a li n e - t o - g r o u n d f a u l t . 9 Inductance between nodes k and m. 31 10 Schematic r e p r e s e n t a t i o n of l o s s l e s s , s h o r t - c i r c u i t e d 32 t r a n s m i s s i o n l i n e . 11 R e l a t i v e amplitude e r r o r of the t r a p e z o i d a l r u l e of 34 i n t e g r a t i o n . 12 Phase e r r o r of the t r a p e z o i d a l r u l e of i n t e g r a t i o n . 35 13 Schematic r e p r e s e n t a t i o n of a three-phase Thevenin e q u i - 40 v a l e n t c i r c u i t . 14 Schematic r e p r e s e n t a t i o n of two connected networks. 45 15 Flow chart of s o l u t i o n w i t h method I . 49 16 Flow chart of the i t e r a t i o n scheme. 55 17 Comparison of f i e l d current w i t h various p r e d i c t i o n 58 techniques. 18 Comparison of s t a t o r current i n phase "a" w i t h v a r i o u s 59 p r e d i c t i o n techniques. 19 Line-to-ground f a u l t at generator t e r m i n a l . 61 20 Simulated current I, i n the u n f a u l t e d phase "b". 62 v i Figure Page 21 System diagram. 63 22 Comparison of s t a t o r currents between s i m u l a t i o n and 64 f i e l d t e s t f o r a three-phase f a u l t . 23 System diagram. 66 . 24 Comparison of the simulated and measured f i e l d c u r r e n t . 66 25 Comparison of the simulated current i n the f a u l t e d phase 67 "a". 26 Comparison of the simulated f i e l d c u r r e n t . 68 27 System diagram. 69 28 Comparison of the simulated current f o r f a u l t e d phase "a". 70 29 Comparison of the simulated three-phase instantaneous 71 power. 30 Comparison of the simulated f i e l d c u rrent. 72 31 Phasor diagram of a synchronous generator f o r balanced, 74 steady s t a t e o p eration. 32 T y p i c a l o p e n - c i r c u i t c h a r a c t e r i s t i c . 81 33 L i n e a r i z a t i o n through the o r i g i n . 85 34 Schematic r e p r e s e n t a t i o n of s a t u r a t i o n e f f e c t s . 88 35 S t r a i g h t - l i n e approximation of the f l u x - c u r r e n t 89 c h a r a c t e r i s t i c . 3.1 Schematic r e p r e s e n t a t i o n of the network. 105 6.1 Flow chart of computer program f o r m a t r i x r e d u c t i o n . 115 7.1 Schematic r e p r e s e n t a t i o n of an unloaded generator.^ 117 8.1 System diagram. 118 8.2 Model of the s h a f t system. 118 8.3 Simulated electromagnetic torque of the generator. 119 8.4 Simulated mechanical speed of the generator r o t o r . 119 8.5 Simulated torque on the s h a f t between the generator 119 r o t o r and the e x c i t e r . v i i LIST OF TABLES Table Page I Comparison of data conversion with different methods. 15 v i i i ACKNOWLEDGEMENT I would like to express my deep appreciation to Dr. H.W. Dommel for his unflagging patience, advice and most importantly faith through a l l phases of this work. Thanks are also due to Dr. A.C. Soudack for his coopera-tion in this project. The financial support from the National Research Council of Canada and the University of British Columbia is gratefully acknowledged. ix 1. INTRODUCTION The importance of general-purpose computer programs for the simulation of electromagnetic transients i n power systems i s constantly increasing. Some of the elements of power systems can now be represented with a high degree of sophistication, e.g., overhead lines with frequency dependent line parameters [1], Some other elements, however, are not yet represented in enough detail, including synchronous generators. Normally, sinusoidal voltage sources E" cos (cot + p ) behind impedances R + jwL'1 6 max y r a J d have been used to represent generators in transient studies. In the derivation of this approximate model, i t i s assumed that rotor fluxes do not change immediately after the disturbance, and that subtransient saliency of the generator can be ignored. This simple model is quite adequate for certain types of studies during the f i r s t cycle or so after the disturbance which initiates the transient phenomena, e.g., switching surge studies, transient recovery voltage studies, and other types of studies involving fast transients [2]. It i s also adequate i f the gen-erator impedance i s only a small part of the total impedance between the generator and the location of the disturbance. Recent interest in electromagnetic transient phenomena which persist over longer time spans makes i t worthwhile to implement more accurate models for synchronous generators into programs for electro-magnetic transients [2-4]. Potential applications include studies of subsynchronous resonance [4-5], dynamic overvoltages, accurate assessment of damping terms in transient s t a b i l i t y studies (due to d.c. offset, harmonics, and asymmetries in short-circuit currents), and other studies. This thesis discusses the major problems of interfacing gen-erator models with an electromagnetic transients program and proposes new 2. solution techniques. F i r s t l y , the choice of an appropriate generator model and the calculation of i t s parameters are discussed. This discuss-ion and a description of numerical problems in the simultaneous solution of the generator equations and the equations of the connected transmiss-ion network provide the necessary background for the introduction of interfacing techniques - around which the major research effort of the thesis was concentrated. The proposed techniques cover in principle a l l the possible approaches to interface problems. Numerical examples are used to test the validity of the proposed techniques. Some additional problems related to the solution of generator equations, e.g., proposed treatment of saturation effects, are described in the f i n a l chapter of this thesis. The contributions of this thesis to power system analysis con-sist of: (a) a c r i t i c a l review of synchronous generator models and selection of a model appropriate for the simulation of electromagnetic transients, (b) a new method for the calculation of synchronous generator para-meters from test data, (c) a new physical interpretation of the discretization error for the trapezoidal rule of integration applied to series inductances, which shows that the resulting difference equations are exact solutions of equivalent lossless stub lines, (d) the development of two alternative interfacing techniques for solving the generator and network equations simultaneously, with one being similar to a technique developed in industry concurrently 2 a . with the research project of this thesis and the other one being a new technique with less restrictions than the f i r s t one, and (e) an analysis and proposed treatment of saturation effects in the synchronous generator. 3. 2. SYNCHRONOUS GENERATOR MODEL 2.1 General Remarks about Physical Device Modelling In general, the derivation of a mathematical model of any phy-si c a l device consists of the following steps [6]: (1) Selection of a model structure based upon observations and phy-s i c a l knowledge; (2) f i t t i n g of parameters of the chosen model to available data; (3) verification and testing of the model; (4) application of the model to i t s given purpose. The basic decisions are made at the f i r s t stage. It i s , for example, necessary to decide whether the physical device can be treated as a linear system. If so, a linear system of equations (differential or algebraic) is used to describe the basic physical phenomena relevant to the device. Therefore, this stage involves some necessary simplifications of the physical reality. At the next stage, relationships between the parameters of the model and the available data have to be established. At this stage, therefore, some additional simplifications may have to be introduced. The last two stages serve as verification of the developed model, and may result in some changes in the model, i f necessary. It should, therefore, be remembered that any mathematical model of a physi-cal device always involves simplifications of physical reality. 2.2 Model of the Electric Part The generator is assumed to be an "ideal synchronous machine" in the sense of Park's definition [7]. The basic assumptions for this ideal generator can be summarized as follows: (1) Saturation effects are neglected. This allows the application of the s u p e r p o s i t i o n p r i n c i p l e , because the model i s then l i n e a r . N e g l e c t i n g the s a t u r a t i o n e f f e c t s i s a common prac-t i c e i n the theory of a l t e r n a t i n g - c u r r e n t machines [8-9]. Techniques f o r i n c l u d i n g n o n l i n e a r e f f e c t s w i l l be discussed l a t e r on. The magnetic c i r c u i t and a l l r o t o r windings are assumed to be symmetrical both w i t h respect to the d i r e c t a x i s , which l i n e s up w i t h the c e n t e r - l i n e through the f i e l d p o l e s , and to the quadrature a x i s 90° behind i t (the recommended p o s i t i o n of the quadrature a x i s l a g g i n g 90° behind the d i r e c t a x i s i s adopted [10 ] ) . A current i n any winding i s assumed to set up a magneto-motive fo r c e s i n u s o i d a l l y d i s t r i b u t e d i n space around the a i r gap. Any magneto-motive f o r c e may be r e s o l v e d i n t o components along the two axes ( d i r e c t and quadrature). The s i n u s o i d a l d i s t r i -b u t i o n does normally imply that only the fundamental component i s considered. In connection w i t h t h i s assumption, i t should be n o t i c e d that the e f f e c t s of harmonics i n the f i e l d d i s t r i -b u t i o n are s m a l l i n a w e l l designed machine [11], [12], I t i s assumed that a magneto-motive forc e a c t i n g along the d i r e c t a x i s produces a s i n u s o i d a l l y d i s t r i b u t e d f l u x wave which also acts along the d i r e c t a x i s . S i m i l a r i l y , a quadrature a x i s magneto-motive forc e produces a s i n u s o i d a l l y d i s t r i b u t e d qua-drature a x i s f l u x . The f a c t o r s r e l a t i n g magneto-motive forc e and f l u x are, however, d i f f e r e n t on the two axes i n a s a l i e n t pole machine [11]. I t i s assumed that the damper bars can be represented as two 5. concentrated hypothetical windings, one in the direct axis (D) and the other in the quadrature axis (Q) [8], Another hypo-thetical winding (g) in the quadrature axis i s normally added for round rotor machines to represent the deep flowing eddy currents. Consequently, the machine consist of seven windings: three a.c. stator windings, one f i e l d winding for direct axis, one hypothetical winding D for direct axis and two hypothetical windings g, Q for the quadrature axis. This "ideal machine" i s schematically shown in Fig. 1. Fig. 1. Schematic representation of a synchronous generator (position of windings is shown in space). A system of seven linear differential equations describes the relationship between the voltages and the currents in the seven windings of the idealized generator. The voltage equations of the generator, in 6. phase c o o r d i n a t e s * , have the f o l l o w i n g form**: where the v e c t o r o f f l u x e s _ i s g i v e n i n g e n e r a l as: _ = [ L ] i = L aa M a b M ac M a f M a D U M ag i a ab ^ b ^ c \f % % M ac V L cc M c f cD cQ M eg i c M a f M . c f L f M f D 0 0 t M a D cD M fD L D 0 0 M « aQ % cQ 0 0 L Q Qg AQ M M eg 0 0 M Q g L g i g (1) (2) The m a t r i x [L] i s always s y m m e t r i c a l , i r r e s p e c t i v e o f r o t o r p o s i t i o n 3. The s e l f and mutual i n d u c t a n c e s o f the armature c o n t a i n even harmonic terms o f r o t o r p o s i t i o n 3 [12], e.g., and L L aa ao L L bb bo L L cc CO C2 ch M = M , + M , eos2B„ + M , , cos4g„ + ... ab abo ab2 3 ahh 3 M = M + M „cos2g + M , cos4g + ... ac aco ac2 2 ac4 2 ^ c " ^ c o + M b c 2 C O s 2 0 l + *W° S 4 61 + (3) (4) * Phase c o o r d i n a t e s r e f e r to a c t u a l c u r r e n t s and v o l t a g e s i n the 3 phases a, b , c o f the s t a t o r and to a c t u a l c u r r e n t s and v o l t a g e s i n the r o t o r w i n d i n g s . ** C a p i t a l l e t t e r s i n square b r a c k e t s [ ] i n d i c a t e m a t r i x q u a n t i t i e s ; s t r a i g h t l i n e s underneath l e t t e r s i n d i c a t e v e c t o r q u a n t i t i e s . 7. where 3i = tot + <S + ir/2 8 2 = Si - 2ir/3 (5) 6 3 = Bj + 2TT/3 The s e l f and mutual inductances of the r o t o r L,., L_, L„, L , f D Q g M£T,, and M„ are constant. The mutual inductances between armature and fD' Qg r o t o r c o n t a i n odd harmonic terms, e.g., M af • M a f l C O s 6 l + M a f 3 C ° s 3 e i + . aD = M a D 1cosB 1 + M a D 3 c o s 3 g 1 + . aQ • M a Q l S l n 6 l + M a Q 3 S l n 3 e i + . M ag " M a g l S ± n e i + M a g 3 s i n 3 3 1 + . (6) S i m i l a r r e l a t i o n s h i p s e x i s t f o r the other two phases b and c [12]. P r a c t i c a l c o n s i d e r a t i o n s a l l o w a re d u c t i o n of the number of parameters appearing i n ( 2 ) - ( 6 ) . For example, a pro p e r l y designed balanced machine i m p l i e s f u l l symmetry of the phases a, b, and c, e.g., (7) L ao = L b o = L CO = L s La2 \ L b 2 = LC2 = L m abo = M aco •^bco = M s ab2 = M „ ac2 = Mbc2 = M 2 Other s i m p l i f y i n g assumptions are made to adjust the complexity of the model to the amount of data which i s u s u a l l y a v a i l a b l e . Some of these can r e s u l t i n n o t i c e a b l e e r r o r s , as demonstrated l a t e r . The f o l l o w i n g l i s t summarizes the most common s i m p l i f i c a t i o n s : (1) A l l harmonic terms of order higher than 2 are neglected [ 8 ] , [12]. This assumption does not seem to cause any n o t i c e a b l e e r r o r s [12]. 8. (2) Mutual inductance M^ (field-to-damper) i s assumed to be equal to the mutual inductances M c- and M (armature-to-field and afl aDl armature-to-damper, respectively), when expressed in p.u. [13], i.e., M = M = M_ = M, afl aDl fD f (8) The same assumption is made for the quadrature axis. These assumptions can sometimes cause significant errors in the simu-lation of rotor quantities [14]. (3) The second harmonic terms in self and mutual inductances of the armature are assumed to be equal, i.e., M2 = L^. This assump-tion simplifies the model significantly, since i t eliminates the coupling among the direct, quadrature, and the zero axes [15]. This simplification is practically always made and seems to be j u s t i f i e d [12], [16], but i t should, nevertheless, be remembered as a possible source of errors. Inclusion of the assumptions mentioned above results in the following, simplified inductance matrix [L]: [L] = L aa Mab M ac M a f M a f Maq ab Lbb ^ c \ f % M ac "be L cc M c f M c f cQ M A cQ M af \ f M cf L f M f 0 0 M af "bf M . cf Mf LD 0 0 M aQ % cQ 0 0 LQ M q M aQ cQ 0 0 M q L g (9) where L = L + L cos28n aa s m 1 L,, = L + L cos2B bb s m 2 L = L + L cos260 cc s m 3 (10) 9. and M , = M + L cos2g„ ab o m 3 M = M + L cos2g„ ac o m 2 M, = M + L cos2g. be o m 1 M . = M.cosg at f 1 \ f = M fcos6 2 M . = M.cosg cf f 3 (ID (12) M n = M s i n g aQ q l % * M q S i n \ M . = M si n g cQ q 3 (13) The r e s i s t a n c e m a t r i x [R] i s simply a diagonal m a t r i x , which has the f o l l o w i n g form f o r a balanced design: [R] = R I R R, R. R g (14) The model of a synchronous generator d e r i v e d above i s b e l i e v e d to be the b e s t , p r e s e n t l y p o s s i b l e compromise between the a v a i l a b l e amount of t e s t data and the d e s i r e d accuracy of s i m u l a t i o n s . R e s u l t s obtained w i t h t h i s model agree q u i t e w e l l w i t h f i e l d t e s t s , which v e r i -f i e s the adequacy of t h i s model [3 ] . The mathematical model of a synchronous generator i n phase-coordinates i s f u l l y d e fined by (1) and (9) - ( 1 4 ) . I t i s , however, common 10. practice in the power industry to describe the generator in a different reference frame, namely in d,q,0-coordinates [7], In this reference frame, a l l inductances defined by (9)-(13) are constant. It should be emphasized, however, that introduction of higher harmonics or unequal second harmonic terms in the armature inductances w i l l result in a time-varying inductance matrix [L] even in d,q,0-coordinates [15]. Before proceeding with transformations to the new reference frame, i t i s useful to rewrite (1) into the following form: v = -[R]i - [L]^= - ~ ( [ L ] ) i (15) The transformation i s defined as follows: and similarly: ip - [ p]'i V = [P]-v ~p — where the subscript "p" denotes Park's d,q,0-quantities. The transformation matrix [P] has the following general form: I | w [P] = 0 where 0 I I L- I J [I] = identity matrix of dimension 4 x 4 : and the matrix [W] is given as [17]: [ W ] = / T c o s B j smt cosB„ sin cosB. sine (16) (17) (18) (19) Park's original transformation matrix and that of many other authors has 2 ^ as a factor, and a negative sign in the second row. The latter i s due 11. to the assumption that the quadrature axis i s leading the direct axis by 90°, rather than lagging behind by 90°. The particular choice of the transformation matrix (19) makes the transformation power invariant and i t s matrix orthogonal, i.e., IP]" 1 = [P] T (20) Application of this transformation to (15) yields: v. = " [P][L]([P] + ^ ( l P r S _ p ) - [P]^([L])[P] ^ - [R]ip -P ,-ld . . d_ dt^p dt x ,d = - [L ]j-± - [R]i - [L ' ] i p dt-p —p P ~V (21) where [L p] = [P][L][P] -1 f M f 0 2 M f 0 0 0 0 0 2 M q 0 # M q ° I M f / I M f ° 0 0 0 L f 0 Mr 0 0 0 0 D 0 0 0 0 \ M 0 0 0 0 M g J (22) and [ L p ] = ^ ^ [ P l ^ t P l t L U P ] h - [ P ] ^ ( [ P ] _ 1 [ L p ] ) = Ui L 0 q -L d 0 0 - / f Mf - / f Mf -J M / 1 M ' 2 q V 2 q 0 0 0 0 0 0 0 0 (23) 0 0 12. 2.3 Calculation of Parameters for the Model of the Electric Part As already explained, any mathematical model of a physical de-vice is an approximation of physical reality. Some of the simplifications made in modelling a synchronous generator were discussed in section 2.2. Here, additional simplifications which are often introduced in the calcu-lation of parameters w i l l be discussed, as well as a technique for avoi-ding them. The functional relationship between measurable parameters, e.g., R , X J } X ' X " , T, ', and T, " i n the direct axis, and the desired ° a d d ' d * do do set of resistances and self and mutual inductances (or reactances) is partly nonlinear [8], [9]. Approximations are normally made to obtain linearized relationships which are easy to solve [11], [18]. These approximations are based upon the knowledge of machine dimensions, which are normally unknown to the system analyst, and were only j u s t i f i e d for hand calculations. There is really no reason any more to introduce them, i f data conversion i s done by a d i g i t a l computer. The model of a synchronous generator was defined in Fig. 1. The inductance matrix [1^] ( i n d,q,0-coordinates) was shown in (22) and the matrix of resistances in (14). As shown in Appendix 1, the measured direct axis machine constants are related to the entries i n the matrices [L p] and [R] in the following way M,2 3 2 Lf " LD ~ 2 M f V = Ld ~ I Mf 2 • <25> L f V M f 1 L f LD 1 / L f LD 2 M f 2 (positive sign of root for T (j 0', negative sign for ") 13. The above equations are nonlinear and their solution by hand is d i f f i c u l t . Introduction of approximations based on the knowledge of machine dimensions leads to a well known, simplified set of equations [11], [19]. The equations (24)-(26) may be solved directly as a system of nonlinear equations by means of Newton's method. Fi r s t , (24)-(26) are rearranged to the following form: G(x) = 0 (27) where G_ = [ g l f g 2, g 3, g^]^ represents the vectors of functional rela-T tionships, x = [R^, L^, L^, R^] represents the unknown machine para-meters, assuming in p.u. is found from: X d = X, + ».Mf (28) Sometimes, the f i e l d resistance R^ in p.u. is given by the manufacturer. In this case, the vector x i s given as: x = [Mf, L f, L D, R^ 1 (29) The following relationship i s obtained by taking the f i r s t term in Taylor's series expansion and equating i t to zero: G(x) = G(x*) + [G'(x*)].(x - x*) = 0 (30) where x represents the approximate solution point of the newest iteration step and x* the approximate solution of the preceding iteration step (or the original guess). From (30) i t follows that x = x* - [G'(x*)] _ 1-G(x*) (31) Eq. (31) i s used iteratively until convergence is reached. The i n i t i a l guesses are found with the approximate linear relationships, as normally used before. This procedure converges quite fast. Typically, only 2 iteration steps are required. 14. A set of equations analogous to (24)-(26) may be obtained for the quadrature axis, which are again solved with Newton's method to obtain the model parameters from the test data. Quite often, the set of data for the quadrature axis is incomplete. In many cases, neither X 1 nor T^o' are given. It is then necessary to reduce the complexity of the model by omitting the g-winding. Equations (24)-(26) are then re-placed by the following set of equations [11] (for details see Appendix 1): M 2 L " = L - |- -9- (32) q 1 Q T " = —2- (33) do RQ Quite often, the manufacturer's data sheets show that T ' ^ T ", but qo qo X^' = Xq» The program would f a i l in this case because the two assump-tions contradict each other (the f i r s t implies the existence of a g-winding, whereas the second implies that there is no g-winding). The following equation: L - M L - L * = M (1 - 9-) = 0 (34) q q q L o would then not have a real solution except for = 0, which implies no g- and Q-winding. This special case can be solved without program modi-fications, however, by setting X^' nearly equal to X^, for example, X ' = 0.99 X (X ' must always be less than X ). Since measurement q q q 3 q accuracy is typically ±5%, this assumption i s acceptable. To il l u s t r a t e the impact of linearization on the calculation of parameters, a 30 MVA machine was considered with the following data [11]: 15. 0.707 p.u. 0.1524 p.u. (calculated from other data) 0.3412 s 314 rad- Is (no g-winding) The r e s u l t s calculated i n three d i f f e r e n t ways are compared i n Table I. Table I. Comparison of data conversion with d i f f e r e n t methods. Conversion method X a q (p.u.) X t (p.u.) X f (p.u.) X D (p.u.) X q (p.u.) S„ (p.u.) R Q (p.u.) a) approximate conversion 1.3879 0.6520 0.0550 0.1795 0.2298 0.1145 0.01766 0.00715 b) exact conversion (with Newton's method) 1.3533 0.6173 0.0897 0.1368 0.1133 0.0698 0.01026 0.00641 c) K i l g o r e ' s method 1.338 0.602 0.105 0.13 0.0678 .0.0515 0.0079 0.0061 While the differences between parameters found from the exact and approximate data conversion methods are not great, the differences i n the simulation r e s u l t s of rotor quantities can be s i g n i f i c a n t . The simulated f i e l d current i ^ i n case of a sing l e line-to-ground f a u l t f o r the generator used i n Table I i s shown i n F i g . 2 f o r approximate and exact data conversion. The i n i t i a l conditions f o r this case are given i n Chapter 3. As pointed out by others [12-14], the res u l t s with approxi-mately calculated parameters have s i g n i f i c a n t l y lower amplitudes of o s c i l -l a t i o n s of rotor q u a n t i t i e s , but the possible improvements with exact parameter conversion has not been recognized. The simulated s t a t o r quantities, V V do rn II do 1.443 p.u. 0.214 p.u. 0.149 p.u. 7.8 s 0.0701 s 0.00064 p.u. X X " q rn II qo Fig. 2. Comparison of the simulated f i e l d currents i ^ . on the other hand, were identical for approximate and exact parameters. This agrees with results published by others [13], [14]. Fig. 3 shows the simulated current 1^ in the faulted phase for the single line-to-ground fault with approximate and exact parameters (indistinguishable from each other). 16.0 r TIME (s) Fig. 3. Identical results for current i c in the faulted phase with approximate and accurate parameter conversion. 2.4 Recent Proposals for Improvements in Parameter Accuracy Additional improvements in the model of a synchronous generator are not possible without introduction of additional test data. New tes-ting procedures require, however, a long time to become accepted as new standards. Model improvements occur, therefore, very slowly, but the data conversion algorithm described above can easily be changed to accommodate additional test data i f and when they become available. The f i r s t improvement i s made i f unequal mutual inductances ^ M^ p k Mp^ are permitted. This was done by a number of researchers. In 1969, Canay [14] explained some of the causes of discrepancies between 18. f i e l d measurements and s i m u l a t i o n r e s u l t s of r o t o r q u a n t i t i e s . He pro-posed an improved e q u i v a l e n t c i r c u i t of the generator, but he re q u i r e d the knowledge of generator dimensions f o r o b t a i n i n g i t s parameters. In 1971, Yu and Moussa [20] suggested another improved model of a generator. They req u i r e d i n t r o d u c t i o n of an e x t r a t e s t (time constant T^ of damper winding) to determine the parameters of t h e i r model. I n 1974, Takeda and Adkins [13] suggested to obt a i n a d d i t i o n a l data from the measurement of the u n i d i r e c t i o n a l f i e l d c u r rent. In the same year, Shackshaft [21] came up w i t h a s i m i l a r , but s l i g h t l y more complicated approach. In a l l cases, the a d d i t i o n a l data can e a s i l y be incorp o r a t e d i n t o the data con-v e r s i o n algorithm, independent of the way i n which i t was obtained. I t i s b e l i e v e d that the data conversion a l g o r i t h m w i l l reduce the d i s c r e -pancies between s i m u l a t i o n and f i e l d t e s t s even more, s i n c e i t avoids the commonly used l i n e a r i z a t i o n of the f u n c t i o n a l r e l a t i o n s h i p s . Some researchers have suggested measuring the parameters d i r e c t l y i n phase coordinates [15], [22]. Such methods w i l l not com-p l e t e l y solve the problems, s i n c e some windings are i n a c c e s s i b l e . They may, however, be very u s e f u l f o r improving the accuracy of the model. A d i f f e r e n t approach has evolved i n the l a s t few years. I t i s based upon parameter e s t i m a t i o n e i t h e r i n the time domain [12] or i n the frequency domain [23], [24], I t i s t h i s author's opinion that the parameter e s t i -mation method w i l l e v e n t u a l l y replace a l l the other approaches. I t may, however, be a long time before t h i s happens. I t was, therefore, b e l i e v e d to be important to develop a simple method making e f f i c i e n t use of the e a s i l y a v a i l a b l e , conventional s et of t e s t data. 19. 2.5 Model of the Mechanical Part of the Generator It is common practice in transient s t a b i l i t y studies to re-present the rotating part of the turbine-generator unit (its shaft, gen-erator and turbine rotors) as one lumped mass. This approach i s , however, unacceptable in some studies of transient performance of the generator, where the rotational vibrations of different parts of the shaft system are important [3], [4]. As in the case of the electrical part of the generator, the complexity of the model depends on the amount of availa-ble data. The actual number of lumped masses may, therefore, vary from case to case. The techniques for modelling mechanical systems of rota-ting lumped masses are relatively well developed. The description of the mechanical part w i l l , therefore, be relatively short. It should be noted, however, that some of the mechanical parameters are very d i f f i c u l t to obtain, not unlike some of the electrical parameters. Fig. 4 shows a typical example of a turbine-generator unit with seven lumped masses, which is based upon an actual case [3], HP IP LP A LPB LPC GEN EXC Fig. 4. Torsional model of a turbine-generator unit. (HP - high pressure turbine, IP - intermediate pressure, LPA, LPB -low pressure units, GEN - generator, EXC - exciter). The dynamic equations of a rotating mechanical system can be derived from Newton's second law [25]. The following general equation can be written for each rotating mass: 20. dt^ i (35) where J = moment of inertia 6 = angle (rotational displacement) ][ T. = sum of a l l torques acting on the rotating mass i To illustrate the use of eq. (35), the three-mass system of Fig. 5 w i l l be used. Assume that r, L L T-z MASS I K 1 2 MASSH K 2 3 MASSm D 2 3 J 2 , D 2 2 J 3 , D 3 3 1 e 2 1% Fig. 5. Schematic representation of a three rotating masses system. i , i+l D. . i i D i , i + 1 J. l T. l shaft stiffness coefficient (between masses i and i+1); viscous self damping coefficient (of the mass i ) ; viscous mutual damping-coefficient (between masses i and i+1); moment of inertia (of the mass i ) ; angle (of the mass i ) ; external torque (applied to mass i ) . 21. Starting with mass 1, the following three equations can be written: for mass 1, d 2e. de, (36) J i ~ T - + D n d r + D i 2 f c - ( 0 i - e 2 ) + K i 2 ( e i - e 2 ) = T l dt for mass 2, d \ . _ d 62 . _ d J 2 — + D 2 2 ^ " + D12dt: ( e2- ei> + D 2 3 — ( 6 2 - 6 3 ) + K12<e2-91> + = T2 dt" 3dt and for mass 3, d2e de. J3^T + D 3 3 ^ + D2 3t- ( e3- 92) + K23( 93- e2) = T3 (37) (38) Equations (36)-(38) assume that the system i s linear, which i s an accept-able simplification for rotational vibrations of small amplitude. A three-mass system has a l l the characteristic features of an n-mass system. Therefore, equations (36)-(38) can be generalized to obtain the following system of equations for any system of n rotating masses: [j] -— e + [ D ] — e + [ K ] e = T dt 2 ~ d t " ~ ~ (39) where [J] = diagonal matrix of moments of inertia; T^ = vector of external torques applied to the system; Q_ = vector of angular displacements; [ D ] = matrix of damping coefficients, which has the following form for the case of Fig. 5, CD] = D N + D 1 2 - D - D 12 12 D 2 2 + D 1 2 + D 2 3 0 - D - D 23 23 D3 3+ D 2 3 (40) 22. [K] = m a t r i x of s t i f f n e s s c o e f f i c i e n t s , which has the f o l l o w i n g form f o r the case of F i g . 5, K 1 2 " K 1 2 . 0 [K] = -K 1 2 K 1 2 + K 2 3 " K 2 3 (41) 0 " K 2 3 K 2 3 In the case of turbine-generator u n i t s , the e x t e r n a l torques are of two types: (1) mechanical i n the t u r b i n e stages; (2) electromagnetic i n the generator and e x c i t e r r o t o r s . The c a l c u l a t i o n of mechanical torques i n t r a n s i e n t s t a b i l i t y s t u d i e s can be simple or very complicated. In the former case, i t i s assumed that mechanical torque or power remains constant a f t e r the disturbance. In the l a t t e r case, the dynamics of the speed governor and a s s o c i a t e d c o n t r o l systems must a l s o be modelled [26], [27]. In electromagnetic t r a n s i e n t s t u d i e s , which i s the s u b j e c t of t h i s t h e s i s , the i s s u e i s l e s s compli-cated due to much shor t e r time spans i n v o l v e d (normally, cases are only simulated up to 1 sec. a f t e r the di s t u r b a n c e ) . I t i s then p o s s i b l e to ass constant mechanical power input and c a l c u l a t e the torque from the f o l l o w -i n g r e l a t i o n s h i p : T = (42) where P = mechanical power p r i o r to disturbance, _ m = angular speed of the mechanical system. I t has been shown that the assumption of constant mechanical torque pro-duces s a t i s f a c t o r y r e s u l t s [ 4 ] , but constant mechanical power seems to be a more reasonable assumption than constant mechanical torque. 23. The electromagnetic torque developed i n the r o t o r of a syn-chronous generator i s equal to the a i r gap torque produced by the r o t a -t i n g electromagnetic f i e l d [9] and may be described by the f o l l o w i n g formulas: (a) i n Park's d,q,0-coordinates [ 9 ] : T e - < V q " Vd } ' f ( 4 3 > (b) i n phase-coordinates [28]: T e = =^ f * a ( 1 b " V + " V + * c ( i a " V ] ( 4 4 ) where n = number of poles of the generator. The torque i n the e x c i t e r ( i f i t i s a d.c. generator d i r e c t l y coupled to the turbine-generator s h a f t ; not modelled i n other cases such as motor-driv e n generators or r e c t i f i e r s ) i s determined by the amount of e l e c t r i c energy produced by the generator and i s given as: T = ( v , - i . + R i b (45) ex 2 »co f f ex f m where co = angular speed of the e x c i t e r ; v^ = voltage at e x c i t e r t e r m i n a l s ; i ^ = e x c i t a t i o n c u r r e n t ; R = armature r e s i s t a n c e of the e x c i t e r , ex The electromagnetic torque c a r r i e s a s i g n opposite to that of the mechani-c a l torque, s i n c e i t represents a load to the mechanical system. 2.6 Conclusions An i d e a l i z e d , l i n e a r model of a synchronous generator, which i s described by r e l a t i v e l y simple equations, has been presented i n s e c t i o n s 0 24. 2.2 and 2.5. The model of the mechanical p a r t was only i n c l u d e d f o r completeness, s i n c e the r e s t of the t h e s i s i s p r i m a r i l y concerned w i t h the e l e c t r i c p a r t of the synchronous generator< This model of a synchronous generator has been shown to be adequate even f o r such complex problems as subsynchronous resonance [ 3 ] , [ 4 ] . I t should be understood, however, that the r o t o r q u a n t i t i e s are not always reproduced a c c u r a t e l y [13], [14]. The problem l i e s i n o b t a i n -i n g enough and s u f f i c i e n t l y accurate data. The numerical c a l c u l a t i o n s , on the othgr hand, can be c a r r i e d out w i t h very high accuracy. The com-p l e x i t y of c a l c u l a t i o n s should, t h e r e f o r e , be r e l a t e d to the accuracy of measurements [29], s i n c e there i s not much sense i n c r e a t i n g a very com-p l e x model f o r ina c c u r a t e data. Lack of r e l i a b l e data o f t e n forces the an a l y s t to use s i m p l i f i c a t i o n s . i n the model. S a t u r a t i o n e f f e c t s have been neglected i n the development of t h i s model. However, as shown l a t e r , i n s e c t i o n 5.5, i t i s p o s s i b l e to in c l u d e them without s a c r i f i c i n g the s i m p l i c i t y of the model. The e l e c t r i c p a r t of a synchronous generator was described i n two systems of coordinates. These two d e s c r i p t i o n s are e q u i v a l e n t f o r t h e o r e t i c a l c o n s i d e r a t i o n s . For numerical s o l u t i o n s , however, one system of coordinates may o f f e r advantages over the other system. This problem w i l l be discussed i n s e c t i o n 3.3. 25. 3. NUMERICAL SOLUTION OF THE GENERATOR EQUATIONS 3.1 Choice of I n t e g r a t i o n Method The dynamic behaviour of a synchronous generator i s described by two sets of d i f f e r e n t i a l equations, one s e t f o r the e l e c t r i c p a r t , and another s et f o r the mechanical p a r t . I t i s important to bear i n mind that no d i g i t a l computer s o l u t i o n of d i f f e r e n t i a l equations can give a continuous h i s t o r y of the t r a n s i e n t phenomena. I t can only give a sequence of "snapshot p i c t u r e s " at d i s c r e t e time i n t e r v a l s At. Such d i s c r e t i z a t i o n causes t r u n c a t i o n e r r o r s , which can lead to numerical i n s t a b i l i t y [30]. The s t e p s i z e At should, therefore,be s m a l l enough t o avoid build-up of t r u n c a t i o n e r r o r s , but not too s m a l l to avoid unnecess-ary computer time. I t i s important to consider the s t r u c t u r a l p r o p e r t i e s of the generator equations i n connection w i t h the choice of the s t e p s i z e At. The system of equations f o r the e l e c t r i c p a r t of the generator and f o r the e l e c t r i c network, to which i t i s connected, i s s t i f f , i . e . the time constants of the system are wid e l y separated [31], In t y p i c a l t r a n s i e n t s t a b i l i t y s t u d i e s , the r a t i o of the l a r g e s t to the s m a l l e s t time constant 3 4 may be i n the order of 10 or 10 [32]. The r a t i o s i n s t u d i e s of e l e c t r o -magnetic t r a n s i e n t s , w i t h which t h i s t h e s i s i s concerned, are s i m i l a r . Even the time constants of the generator equations alone may, i n t h i s rn I case, have a r a t i o i n the order of 10^ or 10"*, e.g., °„ = 1.12*10^ f o r Tdo the generator from s e c t i o n 2.3. In order to avoid numerical i n s t a b i l i t y (due to build-up of t r u n c a t i o n e r r o r s ) , most i n t e g r a t i o n methods, e s p e c i -a l l y those which are e x p l i c i t , r e q u i r e an i n t e g r a t i o n s t e p s i z e At which i s s m a ller than the s m a l l e s t time constant. For i n s t a n c e , fourth-order 26. Runge-Kutta methods require a stepsize which is approximately less than 1/5 of the smallest time constant. Such a small stepsize i s re-quired in spite of the fact that in s t i f f systems the components asso-ciated with the smallest time constants are normally negligible for most of the simulation time span. The overall behaviour of the system, which is of primary interest, i s determined by the largest time constants. The time span of the simulation i s , therefore, determined by the largest time constant. A very small stepsize At i s , therefore, very expensive in simulating s t i f f systems. Round-off errors create additional problems in the numerical solution of di f f e r e n t i a l equations. They may become worse for a solution with a very small stepsize At than for a solution with a larger one. Round-off and truncation error problems are interrelated and are normally considered together as one problem of numerical s t a b i l i t y cf the solution. Any practical method of numerical integration should not only be numerically stable, but also reasonably accurate and e f f i c i e n t . Therefore, the numerical integration method needed in the case of elec-tromagnetic transients should provide a compromise between: (a) numerical s t a b i l i t y ; (b) accuracy; (c) numerical efficiency. The implicit trapezoidal rule of integration seems to be the best compromise for these sometimes contradictory requirements [31], [33]. This method does not suffer from the smallest time constant barrier, i.e., the stepsize At is not controlled by i t . The stepsize At is restricted O mainly by accuracy of the solution, and not by i t s numerical s t a b i l i t y [34]. A fundamental theorem due to Dahlquist [35] states: Theorem: Let a method be called A-stable, i f , when i t is applied to dy the problem -~ = Ay, X < 0, i t is stable for a l l At > 0. Then: (1) No explicit linear multistep method is A-stable; (2) no implicit linear multistep method of order greater than two is A-stable; (3) the most accurate A-stable linear multistep method of order two is the trapezoidal rule: y(t+At)-y(t) = ^|{f(t,y(t)) + f(t+At,y(t+At))} (46) for an equation = f(t,y) (47) The A-stability property was the main reason for the choice of this par-ticular integration method. Some additional, important facts speak in favour of the trape-zoidal rule. F i r s t of a l l , i t i s very simple to program, and does not require past history points except for those of the immediately preceding time step. It i s , therefore, self-starting. It i s also important to note that the trapezoidal rule with a constant stepsize At creates con-stant state transition matrices for linear systems with constant coeffic-ients. This property reduces significantly the amount of calculations involved in the solution process. Finally, i t is worth mentioning that the use of this integration method assures consistency with the Transient Program [34], which uses the same solution method. A number of different solution techniques were suggested in the literature, but none of them seems to have clear advantages over the trapezoidal rule [3], [32], [36-39], Two cases were run to compare the t r a p e z o i d a l r u l e w i t h the fourth-order Runge-Kutta method f o r generator equations. The generator described i n s e c t i o n 2.3 was used i n both cases. F i r s t , a three-phase s h o r t - c i r c u i t was simulated as shown i n F i g . 6. The voltage of the i n -f i n i t e busbar was 2.0/0° p.u., and the i n i t i a l c o n ditions of the genera-t o r were 1.734/-5.2° p.u. s t a t o r current and 3.56 p.u. f i e l d c u r r e n t . The network parameters were R q = 1.0 p.u. F i g . 6. Three phase-to-ground f a u l t at generator t e r m i n a l s . The s i m u l a t i o n r e s u l t s f o r a time step of At = 100 us were p r a c t i c a l l y i d e n t i c a l f o r the fourth-order Runge-Kutta method and the t r a p e z o i d a l r u l e of i n t e g r a t i o n , but some d i f f e r e n c e s were v i s i b l e f o r At = 1 ms. F i g . 7 compares the r e s u l t s f o r the f i e l d current i ^ . As a second t e s t case, a lin e - t o - g r o u n d f a u l t of an unloaded generator was chosen, w i t h the same generator data as f o r the f i r s t case. The fourth-order Runge-Kutta method became n u m e r i c a l l y unstable i n t h i s case, even f o r At = 100 us. The reason may have been the way i n which the e x t e r n a l network was simulated, namely as a very l a r g e r e s i s t a n c e R 29. 0.0 TRAPEZOIDAL 8 RUNGE-KUTTA, At=100ps TRAPEZOIDAL RULE WITH At = 1ms RUNGE-KUTTA WITH At = 1ms 0.0 0-02 0.04 0.06 TIME fs) Fig. 7. Comparison of simulation results for the f i e l d current in case of a three-phase fault. in the two unfaulted phases and zero resistance in the faulted phase. The comparison of the results for two different stepsizes At is shown in Fig. 8. Figs.7 and 8 show that an increase in the stepsize At results in decreased accuracy of the solution with the trapezoidal rule, but not in numerical in s t a b i l i t y . Stepsizes greater than 1 ms were not considered, since the time step for the simulation of r e a l i s t i c cases is dictated by the solution of electromagnetic transients in the external network, where stepsizes are typically in the order of 50 to 100 us. 30. Fig. 8. Comparison of simulation results for the f i e l d current in case of a line-to-ground fault. 3.2 Physical Interpretation of the Trapezoidal Rule of Integration for a Lumped Inductance It i s very important to understand the build-up of discretiza-tion errors when the differential equations of the generator are replaced by difference equations*. As shown earlier, the generator model consists of lumped inductances and resistances. The resistance part does not * The trapezoidal rule of integration applied to v = L^- is identical with replacement of the derivative by a central difference quotient. 31. cause discretization errors, since the functional relationship v = -Ri (48) is solved accurately (except for round-off errors) as a linear algebraic equation. Discretization errors must only be considered for the induc-tance part. To keep the explanation simple, consider a single inductance L between two nodes "k" and "m" as shown in Fig. 9: V m Fig. 9. Inductance between nodes k and m. Application of the trapezoidal rule to yields: i(t+At) = i( t ) + ||{v(t+At) + v(t)} (50) It w i l l now be shown that replacing the differential equation (49) by the difference equation (50) is identical to replacing the lumped inductance L by a short-circuited, lossless transmission line of travel time x = and characteristic impedance Z = This line, which replaces L, has (unavoidable) shunt capacitance C' per unit length which goes to zero as At goes to zero, and a series inductance L' per unit length which, when multiplied with line length, is equal to the value L of the lumped inductance. It is schematically shown in Fig. 10; The equations of the lossless line can be solved exactly with Bergeron's method* [34]. With this method, the following equation can be derived * Method of characteristics in mathematical references. 32. k m H I - - -VS=0 F i g . 10. Schematic r e p r e s e n t a t i o n of a l o s s l e s s , s h o r t - c i r c u i t e d transmission l i n e . f o r a f i c t i c i o u s observer t r a v e l l i n g from t e r m i n a l "1" to "2" i f he leaves t e r m i n a l "1" at t - A t , At, VjCt-At) + Z i ^ t - A t ) = - Z i 2 ( t - ^ ) , At and f o r an observer l e a v i n g t e r m i n a l "2" towards "1" at t j , Z i 2 ( t - ~ ) = V l ( t ) - Z i ^ t ) Summation of (51) and (52) y i e l d s v x ( t - A t ) + Z i ^ t - A t ) + v : ( t ) - Z i ^ t ) = 0 which can be r e w r i t t e n as ±At) = i j C t - A t ) + ICv^ t) + v x ( t - A t ) ) Equation (54) i s i d e n t i c a l w i t h (50), when 2L where Z i s defined as f o l l o w s : Z = Z = At (51) (52) (53) (54) (55) (56) 33. Equation (54) i s an exact s o l u t i o n f o r a l o s s l e s s t r a n s m i s s i o n l i n e [40], Therefore, the a p p l i c a t i o n of the t r a p e z o i d a l r u l e to a lumped induc-tance L i s e q u i v a l e n t to r e p l a c i n g the lumped inductance w i t h a s h o r t -c i r c u i t e d , l o s s l e s s "stub" l i n e w i t h d i s t r i b u t e d inductance L' = L/length and t r a v e l time T = ^ , which i s then solved a c c u r a t e l y . The e r r o r committed by t h i s replacement can be evaluated q u i t e e a s i l y by c a l c u l a t i n g the input impedance of a s h o r t - c i r c u i t e d l o s s l e s s l i n e . The steady-state behaviour of a single-phase t r a n s m i s s i o n l i n e can be described by the f o l l o w i n g general equations [41]: V : = cosh(Y&)V 2 + s i n h ( Y ^ ) - I 2 Z (57) and V 2 I : = sinh(Y^)-2 £ + c o s h ( Y ^ ) ' I 2 (58) w i t h V and I being phasor values. For a s h o r t - c i r c u i t at t e r m i n a l "2" (V 2 = 0 ) , the input impedance seen from "1" can be described as f o l l o w s : Z I N = j ± = Z t a n h ( Y J D (59) For the l o s s l e s s l i n e , the c h a r a c t e r i s t i c impedance i s : z (60) and the propagation constant i s y = /juL*juC' = ju/L'C* (61) This leads to the f o l l o w i n g r e l a t i o n s h i p : Z I N = J t o L , A T ^ ' tan(-~T) (62) The e r r o r i n the frequency domain can now be assessed by c a l c u l a t i n g the r a t i o of the known true value ^2.um-ped ~ ^a)^J t 0 t* i e c o m P u t e d value Z ^ given by (62): 34. H(_) lumped At-cu ,At-coN = c = 5 - C O t ( r - ) ZIN 2 2 (63) Equation (63) is the ratio of the impedance resulting from the application of the trapezoidal rule of integration to the impedance of the lumped inductance. It i s probably more i l l u s t r a t i v e to consider the relative error G(co) = H(_) - 1.0 (64) instead of the absolute error H(co) . The relative amplitude error | G ( I O ) | is plotted in Fig. 11 as a function of " * ^ » t n e P n a se error argG(to) is plotted in Fig. 12 as a function of the same variable. 0 7T_ 2 3jr 2 27T 5TT 2 3JT Fig. 11. Relative amplitude error of the trapezoidal rule of integration. Figs. 11 and 12 give a simple physical explanation of the error for any given stepsize At. Low frequencies are reproduced practically without any error, since lim (|x-cot x - if). = 0 x ->• 0 (65) 35. where x = coAt 2 (66) org Gfoj) degrees 0 TC 2 Tt 3TC 2 2n 57T 2 Tt Y-180 Fig. 12. Phase error of the trapezoidal rule of integration. The error increases with frequency until an absolute blocking (H = 0) is reached for (67) For At = 100 us, a typical stepsize used in studies of electromagnetic transients, the f i r s t blocking frequency is equal to f 1 = 5 kHz. If the frequency is increased beyond f^, then the element is seen by the solu-tion algorithm as i f i t were capacitive up to 2f 1, at which point the element is seen as a short-circuit, reversing to inductive afterwards. The next blocking point i s reached at (68) From there on, the situation repeats i t s e l f periodically. Similar analysis conducted for a lumped capacitance shows that application of the trapezoidal rule to i t s differential equation is equi-valent to replacing the lumped capacitance by an open-ended lossless line. 36. In t h i s case, the r a t i o becomes: H ( w ) = AtT_ t a n ^ > <69) I d e n t i c a l r e s u l t s can be achieved w i t h "the t r a n s f e r f u n c t i o n " approach [42]. The above f a c t s e x p l a i n the sometimes l a r g e l o c a l e r r o r of the s o l u t i o n . The o v e r a l l numerical s t a b i l i t y of the s o l u t i o n , however, i s unaffected [35], [43]. The t r a p e z o i d a l r u l e of i n t e g r a t i o n can thus be regarded as a procedure which changes impedance values d i f f e r e n t l y f o r d i f f e r e n t f r e q u e n c i e s , but solves the system a c c u r a t e l y w i t h these modi-f i e d impedances. The s o l u t i o n i s o b v i o u s l y s t a b l e , but more or l e s s i n a c c u r a t e . 3.3 Choice of Coordinate System In studying electromagnetic t r a n s i e n t s i t i s not c l e a r a p r i o r i which system of coordinates i s more advantageous f o r i n t e g r a t i n g the d i f f e r e n t i a l equations of the e l e c t r i c p a r t . I t i s c l e a r , however, that general-purpose computer programs must i n t e r f a c e the r e s u l t i n g d i f f e r e n c e equations of the generator w i t h those of the network i n phase coordinates. Otherwise the a b i l i t y to simulate any general type of e l e c t r i c network would be l o s t , thus e l i m i n a t i n g the g e n e r a l i t y of the program. The i n t e -g r a t i o n of the generator equations d i r e c t l y i n phase coordinates would seem, t h e r e f o r e , to be the best choice. I t must be r e a l i z e d , however, that the t r a p e z o i d a l r u l e w i l l then produce d i s c r e t i z a t i o n e r r o r s even f o r balanced steady-state c o n d i t i o n s , s i n c e the f l u x and the currents i n the phases change s i n u s o i d a l l y at fundamental frequency. I n t e g r a t i o n i n d^q,0-coordinates, on the other hand, would be exact f o r balanced steady-s t a t e c o n d i t i o n s w i t h the t r a p e z o i d a l r u l e , s i n c e f l u x and currents i n 37. d,q,O-coordinates are constant i n t h i s case. I f the l a t t e r approach i s used, then the r e s u l t i n g d i f f e r e n c e equations must be transformed back to phase coordinates before they are i n t e r f a c e d w i t h the d i f f e r e n c e equa-t i o n s of the network. For general t r a n s i e n t c o n d i t i o n s , the choice i s not at a l l c l e a r . As an example, a s i n g l e l i n e - t o - g r o u n d f a u l t at the generator terminals has phase currents v a r y i n g at 60 Hz, i f the harmonics and dc-o f f s e t are ignored. Currents i n d,q,0-coordinates, on the other hand, w i l l vary at 60 Hz and 120 Hz, w i t h the l a t t e r caused by negative se-quence components. The f i n a l choice was made a f t e r c o n s i d e r i n g the numerical e f f i -ciency. The amount of c a l c u l a t i o n s i n d,q,0-coordinates i s s i g n i f i c a n t l y s m a l l e r , s i n c e the inductance m a t r i x [L ] i s constant i n t h i s case [44]. P P r e l i m i n a r y experiments showed that i n t e g r a t i o n i n d,q,0-coordinates gives very s a t i s f a c t o r y answers, and d,q,0-coordinates were f i n a l l y chosen f o r the work i n t h i s t h e s i s . As shown i n Appendix 2, the t r a p e z o i d a l r u l e of i n t e g r a t i o n i n both systems of coordinates leads to the same form of l i n e a r r e l a t i o n s h i p s between voltages and c u r r e n t s , and only the d i s c r e t i z a t i o n e r r o r s are d i f f e r e n t . The choice of system of coordinates, t h e r e f o r e , does not change the general approach o u t l i n e d i n t h i s t h e s i s . I f space harmonics i n the magnetic f i e l d d i s t r i b u t i o n are taken i n t o account, then i n t e g r a t i o n i n phase coordinates w i l l probably be more advantageous. S e l f and mutual inductances could then be defined d i r e c t l y i n phase coordinates. The inductance m a t r i x [L^] i n d,q,0-coordinates would no longer be constant, i n t h i s case anyhow, thus d i m i n i s h i n g the main advantage of t h i s system of coordinates. 38. 3.4 Multiphase E q u i v a l e n t Networks The a n a l y s i s of connected subnetworks becomes simpler i f each subnetwork has already been reduced to a multiphase e q u i v a l e n t c i r c u i t i n which only terminals at the connection p o i n t s are r e t a i n e d and a l l other terminals are e l i m i n a t e d . Such a s i t u a t i o n e x i s t s when generators are connected to a t r a n s m i s s i o n system. On the generator s i d e , only the s t a t o r windings are d i r e c t l y connected to the transmission network, i . e . , only three p a i r s of i t s terminals must be r e t a i n e d i f the generator equa-t i o n s are to be solved simultaneously w i t h the network equations. This assumes, of course, that the generator equations can i n f a c t be reduced to equations c o n t a i n i n g the r e t a i n e d t e r m i n a l s , only. This i s indeed p o s s i b l e , as explained i n chapter 3.5. A l l other p a i r s of t e r m i n a l s on the r o t o r are then concealed. The s i t u a t i o n i s s i m i l a r i n the t r a n s -m i s s i o n network. Again, only three p a i r s of i t s terminals are connected to each generator, thus making a l l the r e s t of them concealed, provided that the network can be reduced to the r e t a i n e d terminals only. Since the i n t e r f a c e between the electromagnetic t r a n s i e n t s , which solves the network, and the generator model i s performed only through the r e t a i n e d p a i r s of t e r m i n a l s , i t seems appropriate to discuss i n some d e t a i l the i d e a of a multiphase e q u i v a l e n t network. Multiphase e q u i v a l e n t networks have not been used f o r a very long time. Since a number of good references are a v a i l a b l e [45], [46], only a short o u t l i n e w i l l be given here to a i d i n the understanding of i n t e r f a c e techniques. The theory w i l l f i r s t be explained f o r steady s t a t e , w i t h voltages and currents being phasors, and then extended to the s o l u t i o n of electromagnetic t r a n s i e n t s i n s e c t i o n 3.5. This extension becomes s t r a i g h t f o r w a r d w i t h the concept of r e s i s t i v e companion models [47]. 39. The generator or the network can either be described by a set of independent node equations with an admittance matrix, [Y]^ V = I_, or with an impedance matrix, V = [Z]I_ where [Z] = [Y] ^. Retained terminal pairs are generally located across only a few independent node pairs. The node impedance or admittance matrix of the reduced network, which contains only the retained terminal pairs, can be obtained from the f u l l impedance or admittance matrix by elimination of the concealed variables. Consider a general network with N independent node pairs and with R terminal pairs to be retained. Such a network may be described by the following nodal equations with an admittance matrix: Y Y _ _ RR RE "1 Y Y ER — EE _ R^ 2_ (70) or by an inverse relationship with an impedance matrix. ^R ZRR I ZRE -t ZER I ZEE (71) where subscript "R" denotes the retained variables, and subscript "E" denotes the concealed variables, which are to be eliminated. Elimination of the concealed variables results in the following rela-tionships : ,-1 -1 and iR= - [ Y R E H Y E E ] ^ [ Y E R ] ) V R + n^m™] - i ^R= ( [ ZRR ] - r Z R E ] [ Z E E r l [ Z E R ] ) i R + ^ W ^ h (72) (73) Equations (72) and (73) may be interpreted as generalized forms of Norton's or Thevenin's theorem. Equation (73) w i l l be considered closer. It can be rewritten as follows: 40. S e t t i n g I,, = 0. gives the o p e n - c i r c u i t t e r m i n a l voltages as: \ = [ Z R E ] t Z E E ] ~ \ (74) (75) The e q u i v a l e n t impedance [Z ] i s defined as f o l l o w s : = t ZRR ] - ^ R E ^ E E 1 " 1 ^ (76) Equation (74) d e s c r i b e s , t h e r e f o r e , the multiphase Thevenin equivalent c i r c u i t of a network. I t i s a reduced network w i t h concealed terminals e l i m i n a t e d . A n a l y s i s i s s i m p l i f i e d i f the impedance matri x [Z ] behind R the voltage sources i s constant. The Thevenin e q u i v a l e n t c i r c u i t of RU (73) i s only u s e f u l i f the voltages V,, across the concealed node p a i r s are known, and (72) i s only u s e f u l i f the currents 1^ , across the cone cealed node p a i r s are known. A three-phase Thevenin equivalent c i r c u i t i s shown s c h e m a t i c a l l y i n F i g . 13. z R F i g . 13. Schematic r e p r e s e n t a t i o n of a three-phase Thevenin equivalent c i r c u i t . The multiphase e q u i v a l e n t c i r c u i t was derived from nodal equa-t i o n s , but i t can also be derived from other forms of network d e s c r i p t i o n , e.g., from branch equations or mesh equations. The concept of multiphase Thevenin equivalent c i r c u i t s , which was derived f o r steady s t a t e above, can a l s o be used f o r t r a n s i e n t con-d i t i o n s . The e f f i c i e n t c a l c u l a t i o n of a Thevenin eq u i v a l e n t c i r c u i t f o r the transmission network has already been explained i n Appendix 3 f o r t r a n s i e n t c o n d i t i o n s . Therefore, only the Thevenin eq u i v a l e n t c i r c u i t of the generator must s t i l l be derived f o r t r a n s i e n t c o n d i t i o n s , which i s done i n the f o l l o w i n g s e c t i o n . . 3.5 Three-Phase Eq u i v a l e n t C i r c u i t of the Generator A f t e r the a p p l i c a t i o n of the t r a p e z o i d a l r u l e , the nodal v o l t -age equations of a generator i n d,q,0-coordinates have the f o l l o w i n g form ( f o r d e t a i l s see Appendix 2 ) : v ( t ) = [ R C O m p ] i ( t ) + e (t-At) (77) —-p —p —p Equation (77) can be v i s u a l i z e d as v o l t a g e sources e^t-At) behind r e -s i s t a n c e s [ R C O m ^ ] . Such e q u i v a l e n t ' r e s i s t i v e networks, which r e s u l t from the i m p l i c i t i n t e g r a t i o n of d i f f e r e n t i a l equations, are c a l l e d " r e s i s t i v e companion network models" i n network theory [47]. They have been used i n power systems a n a l y s i s f o r more than 10 years [48]. The r e s i s t i v e m a t r i x [ R C O m p ] i s constant i n d,q,0-coordinates, and the v o l tage sources e^(t-At) are c a l c u l a t e d from the known " p a s t - h i s t o r y terms" of the preceding time step t-At. The a b i l i t y to create such " r e s i s t i v e companion network models" i s not l i m i t e d to the t r a p e z o i d a l r u l e of i n t e g r a t i o n only, but works f o r any i m p l i c i t i n t e g r a t i o n as shown i n Appendix 4. Equation (77) represents a system of seven equations. The f i r s t three of them describe the s t a t o r windings, and the r e s t describe the r o t o r windings. Therefore, (77) can be r e w r i t t e n as f o l l o w s 42. ( s u b s c r i p t "p" and s u p e r s c r i p t "comp" dropped to s i m p l i f y n o t a t i o n ) : v (t) = [R ] i (t) + [R ] i (t) + e (t-A t ) (78a) —s ss —s s r — r —s v (t) = [R ] i (t) + [R ] i (t) + e (t-At) (78b) — r r s —s r r — r — r where s u b s c r i p t " s " denotes s t a t o r q u a n t i t i e s , and " r " r o t o r q u a n t i t i e s . E l i m i n a t i o n of r o t o r currents ^ ( t ) ( r o t o r terminals are concealed) r e s u l t s i n the f o l l o w i n g r e l a t i o n s h i p : v ^ t ) = ( [ R s s ] - [ R g r ] [ y 1 [ R „ ] ) i B ( t ) + ' e^(t-At) + [ R s r ] [ R r r ] " 1 . ( e r ( t - A t ) - v ^ t ) ) (79) As mentioned i n s e c t i o n 3.4, the eq u i v a l e n t c i r c u i t of (79) i s r e a l l y only u s e f u l i f concealed v a r i a b l e s are known, which i s v^_(t) here. Since a l l damper windings are permanently s h o r t - c i r c u i t e d , the voltages across these windings are always equal to zero, and the r e f o r e known. Only the f i e l d winding r e q u i r e s s p e c i a l a t t e n t i o n . Depending on the type of study, three approaches can be used: (a) For many st u d i e s t i s so short that the e x c i t e r output does J max not change w i t h i n t h a t time span. The voltage across the f i e l d winding v ^ ( t ) i s then constant and equal to the pre-disturbance value . (b) For st u d i e s over longer time spans, the response of the e x c i -t a t i o n system may have to be taken i n t o account. D i f f e r e n t i a l equations are then used to describe r e l a t i o n s h i p s between t e r m i n a l v o l t a g e , v o l t a g e output v^ of the e x c i t e r , and p o s s i b l y supplementary v a r i a b l e s such as s h a f t speed, a c c e l e r a t i o n , e l e c t r i c power, e t c . I f i m p l i c i t i n t e g r a t i o n i s a p p l i e d to these d i f f e r e n t i a l equations ( l i n e a r or l i n e a r i z e d over one 43. time s t e p ) , then l i n e a r a l g e b r a i c r e l a t i o n s h i p s among these q u a n t i t i e s are obtained. With a t y p i c a l time constant of 30 ms f o r a f a s t e x c i t a t i o n system w i t h r e c t i f i e r s , i t should be p e r m i s s i b l e to c a l c u l a t e the e x c i t a t i o n system output at time t from the known input values at time t - A t , and v ^ ( t ) would then again be known i n the s o l u t i o n of the generator equations. This approach was used s u c c e s s f u l l y f o r p r a c t i c a l cases [49]. Standard IEEE e x c i t a t i o n system models define the t e r m i n a l v o l t a g e as an RMS-value. Therefore, the problem a r i s e s i n t r a n s i e n t s t u d i e s how to d e f i n e RMS-values from instantaneous values under t r a n s i e n t c o n d i t i o n s . This could, f o r example, be done by i n c l u d i n g a model f o r the t r a n s i e n t behaviour of the transducer. This i s s u e may r e q u i r e f u r t h e r research, (c) I f the time delay of one time step introduced i n method (b) above i s unacceptable, then i t becomes necessary to r e t a i n the f i e l d winding i n the e q u i v a l e n t c i r c u i t , which leads to a f o u r -phase e q u i v a l e n t c i r c u i t of the generator. The f o u r t h equa-t i o n must then be i n t e r f a c e d w i t h the equations which describes the e x c i t a t i o n system dynamics. This w i l l not a f f e c t the i n t e r f a c e procedures described i n t h i s t h e s i s ( f o r d e t a i l s see Appendix 5 ) . Equation (79) can now be r e w r i t t e n as f o l l o w s : v ( t ) = [ R r e d ] i (t) + e r e d ( t - A t ) —s ss —s —s (80) w i t h [ R r 6 d ] = ss a l l A 1 2 0 A 2 1 A 2 2 0 0 0 A 3 3 (81) 44. where a l l nonzero elements are constant i f At and co are constant ( f o r red d e t a i l s see Appendix 5 ) , and where e (t-At) i s known at time t from s known past h i s t o r y at t-At and v ^ ( t ) . Equation (80), a f t e r transformation to phase-coordinates, describes the three-phase e q u i v a l e n t c i r c u i t of the generator i n those coordinates. The r e s u l t i n g r e s i s t i v e companion m a t r i x [R^™^] i s time dependent: b n ( t ) b 1 2 ( t ) b 1 3 ( t ) [ R c o m P ] = b 2 1 < t ) b 2 2 ( t ) (82) b 3 1 ( t ) b 3 2 ( t ) b 3 3 ( t ) C a l c u l a t i o n of the m a t r i x [R ] according to (79) i s numeri-ss c a l l y very i n e f f i c i e n t and i s b e t t e r done i n p r a c t i c e w i t h a Gauss-Jordan e l i m i n a t i o n scheme [50], as b r i e f l y explained i n Appendix 6. 4. INTERFACING THE GENERATOR MODELS WITH THE TRANSIENTS PROGRAM 4.1 Problem Formulation The generators and the network to which they are connected, can in principle be. solved as one system of equations. However, i t is then not easy to write general-purpose computer programs which can handle any network configuration [51]. It i s , therefore, necessary to devise interfacing procedures which preserve the generality of the net-work representation in the Transients Program. When elect r i c networks are connected together, then certain boundary conditions must be satisfied for voltages and currents at the connection points. The situation for two three-phase networks i s shown in Fig. 14. The conditions, which must be satisfied at any time t, are based upon Kirchhoff's laws, in this case: V l = V V2 = V V3 = \ a n d *•! + \ = °» 12 + i5 = °» "^3 + 16 = °' {SUBNETWORK I '4 >2 '5 SUBNETWORK '3 >6 n Vl V 5 V 6 / / / ) / / Fig. 14. Schematic representation of two connected networks. 46 These conditions must also be satisfied when two computer pro-grams are interfaced, when each of them describes the behaviour of one subnetwork only. Subnetwork I could be the generator, and subnetwork II the transmission system in a particular case. The generator i s connected to the transmission network through three pairs of terminals, i.e., a l l the other pairs of terminals of the transmission system are, from the generator's point of view, concealed. Similarly, looking from the trans-mission system into the generator, only the three stator pairs of gen-erator terminals are visible, and the rotor windings are concealed. There are, therefore, two possible ways of interfacing a generator pro-gram with a network transients program. The f i r s t one is based upon the calculation of a three-phase Thevenin equivalent circuit of the trans-mission network, as seen from the generator's stator terminals, and solving i t together with the f u l l set of generator equations or with the reduced set of generator equations, in which only stator terminals are retained. The f i n a l solution in the transmission network is obtained by superimposing the voltage changes which results from the generator cur-rents on the solution obtained without generator currents [2]. The second approach is based upon the development of a three-phase Thevenin equivalent circuit for the generator in the form of a voltage source red behind a time-invariant, symmetrical matrix [Rp^ ]• Th e complete solu-tion i s then obtained by solving the transmission network with the red generator treated as voltage sources behind ] i n o n e step. This approach results in a significantly simpler interface code. 47. 4.2 Method I - Interface by Means of a Thevenin Equivalent Circuit of the Transmission Network The Transients Program was modified to produce a three-phase Thevenin equivalent circuit of the transmission network seen from the generator terminals (in i t s original form i t could only produce single-phase equivalent circuits). The equivalent circuit can be described by the following equation: VO = [ ^ e r m ± n a ± ] ^ t ) + % Q ( t ) (83) where the subscript "N" denotes network quantities, and the subscript " 0 " denotes open-circuit quantities. The.3 x 3 matrix r R ^ e r m ^ n a ^ ] , i s precomputed before entering the time step loop, and must only be recomputed when the network configuration changes due to switching actions or when the program moves into a new segment in the piecewise linear representation of nonlinear elements [45] A practical way of calculating this matrix, as done in the Transients Program, is briefly described in Appendix 6. The voltage vector y^nCt) contains the three-phase voltages of the Transients Program solution without the generator. The generator is represented by i t s three-phase equivalent c i r c u i t : v ^ C t ) = i R j ^ V t ) + % h ( t " A t ) ( 8 4 ) where the subscript "ph" denotes generator quantities in phase coordin-ates, and where V t ) B ^ ( t ) a n d V 0 " " ^ 0 0 ( 8 5 ) With (85), equations (83) and (84) can be solved for the unknown voltages 48. and c u r r e n t s . A f t e r s u b t r a c t i n g (84) from (83), the f o l l o w i n g r e l a t i o n -s h i p i s obtained: The f i n a l network s o l u t i o n i s found by superimposing the v o l t a g e changes Av(t) - [ l ^ e t W O r k ] - V t > C87) on the previous s o l u t i o n without the generator. Ay(t) i s the v e c t o r of a l l v oltages on nodes without voltage sources i n the t r a n s m i s s i o n net-work, and [ j ^ e t w o r k ] ^ g a p r e C 0 m p u t e d n x 3 m a t r i x from which [ R ^ e r m ^ " n a ^ ] was e x t r a c t e d f o r (84). The f i n a l s o l u t i o n of the generator i s found by s o l v i n g f o r the concealed v a r i a b l e s which were e l i m i n a t e d i n reducing the equations to a three-phase Thevenin eq u i v a l e n t c i r c u i t . A f t e r transformation of the s t a t o r currents to d,q,0-coordinates, the r o t o r currents are found as f o l l o w s : i (t) = -[R ] _ 1[R ] i (t) - [R. ] _ 1 ( v (t) + e ( t - A t ) (88) — r r r 1 r s —s r r J v — r — r -1 -1 where'the matrices [R ] and [R 1 [R ] were found as by-products of r r r r r s the r e d u c t i o n process, as explained i n Appendix 6. A flow chart f o r t h i s s o l u t i o n a l g o r i t h m , w i t h the mechanical p a r t of the turbine-generator i n c l u d e d , i s shown i n F i g . 15. The p r e d i c t i o n p a r t of the s o l u t i o n a l g o r i t h m i s only needed i f the mechanical p a r t i s i n c l u d e d f o r modelling the r e l a t i v e r o t o r o s c i l l a t i o n s around synchronous speed. Since the i n t e r f a c i n g i s done i n phase coordinates, i t i s necessary to know both the angular p o s i t i o n 0 and the speed u of the r o t o r at the new time step i n order to c a l c u l a t e m the m a t r i x [R C° m p] and the v e c t o r e , ( t - A t ) . A s i m i l a r s i t u a t i o n e x i s t s pn —ph START : t=0 t = t+At 49. TRANSIENTS ~»?•:>.,\y. SOLUTION ' J ITKOUT GENERATORS CALCULATE H I S T . TERMS I PREDICT ROTOR ANGLES 6 AND SPEEDS 1% CALCULATE CURRENTS, ROTOR ANGLES 6 AND SPEEDS ium WITH TERMINAL INFORMATION FROM TRANSIENTS PROGRAM YES SUPERIMPOSE VOLTAGES PRODUCED BY GENERATOR TO THOSE OBTAINED BEFORE NO | STOP | •J Z o to pes O H te Fig. 15. Flow chart of solution with method I. 50. when the i n t e r f a c i n g i s done i n d,q,0-coordinates. I t i s then necessary to know 6 and OJ i n order to c a l c u l a t e the m a t r i x [ R C O m p ] and the three-m p phase Thevenin equivalent c i r c u i t of the t r a n s m i s s i o n system i n Park's coordinates. A f t e r the generator currents have been c a l c u l a t e d , i t i s p o s s i b l e to c a l c u l a t e the electromagnetic torque T£^. With t h i s v a lue, the equations of the mechanical p a r t can be solved to get updated value of the speed t o ^ . I f these values d i f f e r too much from the p r e d i c t e d v a l u e s , then the s o l u t i o n i s repeated u n t i l the d i f f e r e n c e s are n e g l i g i -b l e . This s o l u t i o n a l g o r i t h m performed s a t i s f a c t o r i l y [4]. C o r r e c t i o n s are not much of a problem i n t h i s s o l u t i o n a l g o r i t h m , anyhow. 4.3 L i m i t a t i o n s of Method I The s o l u t i o n method described i n chapter 4.2 i s q u i t e s t r a i g h t -forward and n u m e r i c a l l y s t a b l e i f only one generator i s connected to the t r a n s m i s s i o n network. I f there are more generators connected to the network, as i s u s u a l l y the case, then the method works without f u r t h e r m o d i f i c a t i o n s only i f the generators are separated through d i s t r i b u t e d -parameter l i n e s , i . e . , i f three-phase e q u i v a l e n t c i r c u i t s of the net-work e x i s t independently f o r each generator because distributed-parameter l i n e s disconnect the network [ 2 ] , [44]. The Transients Program checks t h i s c o n d i t i o n a u t o m a t i c a l l y when i t c a l c u l a t e s three-phase e q u i v a l e n t c i r c u i t s at generator t e r m i n a l s . The program could be changed to c a l c u -l a t e 6-phase equ i v a l e n t c i r c u i t s i f two generators are connected to the network without separation through distributed-parameter l i n e s , or 9-phase eq u i v a l e n t c i r c u i t s f o r three generators, e t c . For the general case of m generators, however, the programming would become.too complicated and the execution time would probably become too b i g . 51. A number of alternative methods are possible to overcome this limitation. The generators could be separated by short distributed-parameter transmission lines. This may, however, require a significant reduction of the stepsize At in order to keep the length of such a r t i -f i c i a l lines short, since At must be less than the travel time [ 4 8 ] , A slightly different approach has been used successfully by Southern California Edison Company, whereby transformer leakage inductances of step-up transformers are approximated as stub-lines [ 4 4 ] . If a group of generators not separated by transmission lines feeds into the same bus-bar, then the possibility exists of creating one equivalent ci r c u i t for this group of generators, which would be relatively easy. This method can be extended to groups of generators which feed through unit trans-formers into the same busbar. In this case, the generator equation would have to be expanded to include the transformer equations as well. The limitations mentioned above do not significantly degrade the practical value of method I. They simply imply that certain pre-cautions are needed when special cases are simulated. 4 . 4 Method II - Interface with a Modified Thevenin Equivalent Circuit of the Generator It i s common practice in the power industry to represent gener-ators by sinusoidal voltage sources of the form E" cos (o)t+p) behind max " subtransient impedances R Q + jooL^| in quasi-steady-state fault studies and in some types of transient studies. In the derivation of this model i t is assumed that flux changes in the rotor windings immediately after the distrubance can be ignored, and that subtransient saliency can also be ignored. This simple model i s quite adequate for the f i r s t cycle or so 52. a f t e r the disturbance which i n i t i a t e s the t r a n s i e n t phenomena. Good r e s u l t s have been obtained w i t h t h i s model i n s w i t c h i n g surge s t u d i e s , t r a n s i e n t recovery v o l t a g e s t u d i e s , and other types of s t u d i e s i n v o l v i n g f a s t t r a n s i e n t s . The r e l a t i v e l y accurate r e s u l t s obtained w i t h t h i s simple model motivated the research e f f o r t described i n t h i s chapter. The i d e a was to f i n d a way to account f o r the changes i n f l u x e s and to i n c l u d e the s u b t r a n s i e n t s a l i e n c y without l o s i n g the s i m p l i c i t y of the model. Before proceeding w i t h the d i s c u s s i o n , i t may be u s e f u l to r e c a l l some of the r e s u l t s of s e c t i o n 3.5, where the r e t a i n e d s t a t o r v a r i a b l e s of the generator were described i n d,q,0-coordinates by the equation v ( t ) = [ R ^ d ] i ( t ) + e ! e d ( t - A t ) (89) —s ss —s . —s red and where the m a t r i x [ R o o ] was given as f o l l o w s : a l l a12 0 ss a22 0 (90) 0 0 a33 As mentioned i n s e c t i o n 3.5, the m a t r i x [ R ^ " ] becomes time-s s dependent when i t i s transformed d i r e c t l y to phase-coordinates. While nodal network s o l u t i o n s , such as i n [34], can i n p r i n c i p l e be m o d i f i e d to accept time-dependent r e s i s t a n c e s , program execution would be slowed down s i g n i f i c a n t l y i f the network conductance m a t r i x had to be r e t r i a n -g u l a r i z e d i n each time step. A l s o , the conductance matr i x becomes asym-m e t r i c i n t h i s case, which means that the upper as w e l l as the lower t r i a n g u l a r matrices would have to be s t o r e d . This would increase s t o r -age requirements as w e l l as computation time compared w i t h the present 53. method based upon symmetry where only the upper t r i a n g u l a r m a t r i x i s s t o r e d . A number of ways were t r i e d out to approximate (89) i n such a way that the r e s i s t a n c e m a t r i x becomes constant and symmetric i n phase coordinates. The f o l l o w i n g scheme seemed to work b e s t : red (1) S p l i t the m a t r i x [R ] i n t o the sum of two terms, s s [ R r e d ] = [ R r e d 1 + [ R r e d ] ss const var (91) red where the m a t r i x [R ] i s given as const ° [ R r e d J -const r- a +a 11 22 2 0 0 a +a 11 22 2 0 33 (92) w i t h c o e f f i c i e n t s a ^ as defined i n Appendix 5. red (2) Transform the m a t r i x [R 1 to phase coordinates. The const 27 fid r e s u l t i s a constant symmetric m a t r i x [ R ^ ] of the f o l l o w i n g form: ,red. b a b b b (93) The elements of the m a t r i x [R J' C U] are normally much sm a l l e r than those var J of [ R ^ ; s t ] , and t h e i r i n f l u e n c e i n (89) can be accounted f o r by m u l t i p l y -i n g them w i t h the p r e d i c t e d values of s t a t o r c u r rents r a t h e r than w i t h a c t u a l , y e t unknown, values and adding these terms to the v o l t a g e sources. Therefore, the f o l l o w i n g r e l a t i o n s h i p i n d,q,0-coordinates i s obtained: v ( t) = [ R r e d J i ( t ) + { e r e d ( t - A t ) + [ R r e d ] i p r e d ( t ) } —s const —s —s var —s v (94) 54. where i ^ r e d ( t ) i s the v e c t o r of the p r e d i c t e d d,q,0-stator currents at time t . Transformation of (94) to phase-coordinates y i e l d s the d e s i r e d Thevenin e q u i v a l e n t c i r c u i t of the generator: ' _p h<t) = [ * p f ] i C t ) + e j ; d ( t ^ t ) (95) red Since the m a t r i x [Rp^ 1 i s symmetric and time-independent, the genera-t o r can be represented i n the Transients Program simply as a s e t of voltage sources e ^ ^ d ( t - A t ) behind r e s i s t a n c e s [ R p j ^ ] . i n v e r s e of red [Rpk ] enters i n t o the nodal conductance m a t r i x of the t r a n s m i s s i o n network as any other r e s i s t i v e branch, as described i n [48], A f t e r the network s o l u t i o n has been found i n the new time step at time t , the s t a t o r q u a n t i t i e s w i l l be known, and the r o t o r currents can then be found from (88). Therefore, the s o l u t i o n of the transmission network together w i t h the generators proceeds at each time step as f o l l o w s ; (1) Solve the network equations together w i t h the Thevenin equiva-l e n t c i r c u i t s of the generators reduced to the s t a t o r windings as given by (95); (2) Solve equation (88) of r o t o r windings, using the s t a t o r currr r e n t s found i n the previous step. There i s a need to p r e d i c t the angular p o s i t i o n 9 and the speed to of the r o t o r , j u s t as i n method I . The problem of i t e r a t i v e c o r r e c -t i o n s i s , more complex i n t h i s case, however, since i t e r a t i o n s w i t h the complete network s o l u t i o n would be q u i t e c o s t l y . F o r t u n a t e l y , i t has already been shown that the angular p o s i t i o n 6 can be p r e d i c t e d a c c u r a t e l y enough without need f o r c o r r e c t i o n s [ 4 ] , [32]. A reasonably accurate p r e d i c t i o n f o r the angular speed _ can probably be obtained from the 55. i n t e g r a t i o n of the'mechanical equations, w i t h l i n e a r e x t r a p o l a t i o n of the e l e c t r i c a l torque because of the r e l a t i v e l y s m a l l s t e p s i z e At used i n i n t e g r a t i o n of the e l e c t r i c p a r t i n comparison w i t h the r e l a t i v e l y b i g time constants of the mechanical p a r t . Furthermore, the terms con-t a i n i n g the angular speed co are only a s m a l l p a r t of e , (t- A t ) i n (95) m —ph I t i s , t h e r e f o r e , reasonable to expect 1 that c o m can be p r e d i c t e d w i t h s u f f i c i e n t accuracy. I t i s al s o p o s s i b l e to introduce some s o r t of i t e r a t i o n loop to improve the s o l u t i o n such as the one shown i n F i g . 16. p r e d i c t 6 and co m No Solve the network together w i t h the generators C a l c u l a t e r o t o r currents and e l e c t r o m a g n e t i c a l torques C a l c u l a t e c o m and 9 w i t h the new data F i g . 16. Flow chart of the i t e r a t i o n scheme. 56. No s e r i o u s problems are expected. Experience gained by the General E l e c t r i c Company seems to prove t h i s statement [ 3 ] . A computer program developed by t h i s company and the r e s u l t s obtained w i t h i t seem to i n d i c a t e that no s i g n i f i c a n t e r r o r s are introduced when using the pre-d i c t e d values of the angular speed to . m 4.5 Remarks about Method I I As already mentioned, the elements of the r e s i s t a n c e m a t r i x red red [R ] are normally much sm a l l e r than the elements of the m a t r i x [R 1. var 3 const This i s true of t y p i c a l generators where s u b t r a n s i e n t s a l i e n c y i s very s m a l l . I f d i f f e r s s i g n i f i c a n t l y from X'"j, the elements of the m a t r i x [ R v a r ] may become l a r g e . This would increase the r e l a t i v e weight of the i n a c c u r a c i e s r e s u l t i n g from the p r e d i c t i o n of the s t a t o r c u r r e n t s , which i n t u rn could r e s u l t i n i n a c c u r a t e s o l u t i o n s . Although such a case seems to be h i g h l y improbable f o r any p r a c t i c a l generator equipped w i t h damper windings, i t should, n e v e r t h e l e s s , be mentioned as a p o s s i b l e problem. Experience has shown t h a t f o r s t e p s i z e s At i n the order of -4 10 s, method I I works remarkably w e l l . I t does not s u f f e r from the l i m i t a t i o n s t y p i c a l of method I , s i n c e the generators are simply modelled as voltage sources behind r e s i s t a n c e s . In the network s o l u t i o n i t i s p o s s i b l e , t h e r e f o r e , to have any number of generators connected to the network, e i t h e r at the same busbar or at d i f f e r e n t busbars, without l o s s of g e n e r a l i t y . The p r e d i c t i o n of the d,q,Q-?stator currents i ^ r e d ( t ) does s i n f l u e n c e the accuracy, of course, and can be performed i n a number of ways: (1) Assume that the voltages at the generator terminals are cons-tant over the next time step, and use the new currents found from (89) as p r e d i c t e d values. (2) Use s t r a i g h t - l i n e e x t r a p o l a t i o n of the c u r r e n t s . Information f o r two preceding time steps at t-2At and t-At must then be s t o r e d . (3) Use p a r a b o l i c e x t r a p o l a t i o n of the c u r r e n t s . Information f o r three preceding time steps must then be s t o r e d . (4) Use any combination of the three previous methods, e.g., s t r a i g h t - l i n e e x t r a p o l a t i o n of the voltages combined w i t h approach 1. In a l l the t e s t s conducted f o r t h i s t h e s i s , there were no v i s i b l e d i f f e r ences between the r e s u l t s obtained w i t h d i f f e r e n t p r e d i c t i o n methods as long as some type of p r e d i c t i o n of the s t a t o r currents was used. Simply s e t t i n g i | \ r e d ( t ) = 0_ i s too i n a c c u r a t e . F i g . 17 compares the f i e l d cur-rent from the example of s e c t i o n 2.3, c a l c u l a t e d i n three d i f f e r e n t ways (a) The p r e d i c t e d currents are obtained by assuming that the v o l -tages at the generator terminals are constant over the next time step, w i t h s t e p s i z e At = 100 us. (b) The currents from the previous time step are used as p r e d i c t e d v a l u e s , i . e . , i p r e d ( t ) = i ( t - A t ) , (At = 100 ys) . s s (c) The i n f l u e n c e of the term [R ] i ( t ) i s neglected, i . e . , var —s i p r e d ( t ) = 0 (At = 100 y s ) . This amounts to n e g l e c t i n g the —s sub t r a n s i e n t s a l i e n c y and speed terms i n the generator model. The d i f f e r e n c e s are l a r g e r f o r the s t a t o r c u r r e n t s , as shown i n F i g . 18 f o r one of the u n f a u l t e d phases (phase a) . ' Using current values from the previous time step t-At can only be j u s t i -f i e d f o r a s m a l l s t e p s i z e At. With i n c r e a s i n g At, the e r r o r introduced by t h i s simple p r e d i c t i o n may become i n t o l e r a b l e , and p o s s i b l y l e a d to numerical i n s t a b i l i t y . 58. 14.0 r 0.0\ , , , , , 0.02 0.04 0.06 TIME fs) Fig. 17. Comparison of f i e l d current with various prediction techniques. The low sensitivity to the accuracy of prediction of the stator currents may serve as an indication of the behaviour of the solution for inaccurately predicted angular speed _ . It is possible to interpret the error in prediction of w as an additional error in the predicted stator currents. It i s , therefore, believed that no serious problems w i l l arise with the introduction of the mechanical part of the generator. red The elements of the matrix 1 c a n ^ e P r e c a l c u l a t e ( i before entering the time step loop. Simple matrix multiplication (transforma-tion from one coordinate system to another) and basic algebra show that: 59. Fig. 18. Comparison of stator current in phase "a" with various prediction techniques. a + a + a a = - i i ? 2 33 ( 9 6 ) and 2 a 3 3 " a i l " a22 b M ^1 22 ( 9 ? ) const where a and b are elements of the matrix [R , ] defined in (93). These Ph elements remain constant as long as the stepsize At does not change and as long as nonlinear saturation effects are ignored. It i s worth mentioning that neither method I nor method II are tied to the trapezoidal rule of integration. The a b i l i t y to create 60. resistive companion network models, and to develop reduced three-phase equivalent circuits from them, i s a general property of any implicit integration scheme, as shown in Appendix 4. The trapezoidal rule of integration, however, makes the numerical procedure simple. This fact, and the s t a b i l i t y properties of this particular integration as discussed in section 3.1 f u l l y justify i t s use. 4.6 Numerical Examples It i s very easy to set up hypothetical generator and network cases for testing the simulation methods of this thesis. While this is satisfactory for comparing various approaches, i t does not answer the question of how closely such simulations agree with f i e l d tests. Every effort was, therefore, made to duplicate published f i e l d test results, even though not too many such cases could be found. In this connection, i t should be remembered that the generator model alone does basically not need verification, since this was already done by others [ 4 ] , [11-14]. The test examples are, therefore, mainly used to verify the various numerical procedures and interfacing techniques. This was done in stages: (1) Preliminary Tests In the preliminary testing of the methods, the interface with the general-purpose Transients Program was replaced by a simple three-phase Thevenin equivalent circuit (voltage source behind external resis-tance and inductance). This simplification made i t easier to test the simulation techniques before interfacing the generator subroutine with the large Transients Program. The program performance was checked for a number of examples. A simulation program using a fourth-order Runge- , Kutta integration routine was also run in parallel as a check on the 61. s o l u t i o n w i t h the t r a p e z o i d a l r u l e . Some of these r e s u l t s were already presented i n s e c t i o n 3.1 ( F i g . 7 and F i g . 8). (2) Example 1 f o r Test i n g Method I The r e s u l t s of the p r e l i m i n a r y t e s t s were compared w i t h the r e s u l t s obtained w i t h method I . This provided a check on both the c a l c u -l a t i o n of the three-phase Thevenin equivalent of the transmission network and on method I i t s e l f . Only one example w i l l be shown, s i n c e there were no v i s i b l e d i f f e r e n c e s . The generator described i n s e c t i o n 2.3 was used w i t h a line-to-ground f a u l t a p p l i e d to one of i t s t e r m i n a l s . The system was simulated as shown i n F i g . 19. The network parameters were R g = 1.0 p R + R = R , and R =0.01 p.u. The voltage of the i n f i n i t e busbar ei ez e el •> was 2.0/0° p.u., and the i n i t i a l c o n d i t i o n s of the generator were 1.7341-5.2° p.u. s t a t o r current and 3.56 p.u. f i e l d c u rrent. Both the s t a t o r currents and the r o t o r currents were i d e n t i c a l f o r the two d i f -f e r e n t s o l u t i o n s (without and w i t h the c a l c u l a t i o n of the Thevenin equi-v a l e n t c i r c u i t ) . The s t a t o r current i n the unfau l t e d phase "b" i s shown i n F i g . 20. INFINITE BUSBAR F i g . 19. Line-to-ground f a u l t at generator t e r m i n a l . 62. F i g . 20. Simulated current 1^ i n the unfa u l t e d phase "b". The simulated f i e l d current was already shown i n F i g . 2 and the f a u l t current I was already shown i n F i g . 3. (3) Example 2 f o r Method I This i s a t e s t case w i t h an e x t e r n a l network f o r which f i e l d t e s t r e s u l t s were a v a i l a b l e , and which could no longer be solved w i t h the program used f o r the p r e l i m i n a r y t e s t s . The generator had the f o l l o w -i n g data (based on 150 MVA and 13.8 kV, RMS, l i n e - t o - l i n e ) [53]: X. V V X., V 1.85 p.u. 0.2575 p.u. 0.18 p.u. 0.175 p.u. 0.85 s 0.385 s X X " q r n II q x 1.76 p.u. 0.29 p.u. (assumed according to [53]) 0.04 s (assumed according to [53]) 0.198 p.u. (no g-winding) 63. The data conversion to [R] and [L p] was done with the formulas published in [11]. The generator was connected to a transmission system, as shown in Fig. 21. 13.8/130-8 k VA . SINGL E PHA SE a CLOSE AT t = o CLOSE AT -i L l h = 10ms TTT Fig. 21. System diagram. The system was i n i t i a l l y unloaded and the voltage at i t s ter-minals was 13.8 kV 7-10° kV (RMS, line-to-line). A three-phase fault on the high side of the delta/wye connected step-up transformer was studied, with the closing sequence as shown in Fig. 21. The simulated stator currents are compared with the f i e l d test results, as given in [54], in Fig. 22. The f i e l d current was unavailable for comparison. The differences in the i n i t i a l values of the currents result from uncertainty in the generator data, e.g., the values for X^" and T^" were assumed rather than measured and may be unrealistic, and lack of sufficient information about the i n i t i a l conditions, e.g., about the instant at which the fault was applied. Because of the latter reason, the i n i t i a l conditions had to be varied until the results came reasonably 64. (c) Fig. 22. Comparison of stator currents between simulation and f i e l d test for a three-phase fault, (a) = phase "a", (b) = phase "b", (c) = phase "c", — f i e l d test, — simulation. 65. c l o s e to those of the f i e l d t e s t s . I t should a l s o be remembered t h a t s a t u r a t i o n e f f e c t s were ignored. In view of a l l these l i m i t a t i o n s the agreement i s reasonably good. (4) Example 3 f o r Tes t i n g Method I I t was found by a number of researchers t h a t the c o r r e c t r e -production of s t a t o r currents i s not much of a problem [13], [14]. I t i s o f t e n d i f f i c u l t , however, to reproduce r o t o r q u a n t i t i e s c o r r e c t l y . To i l l u s t r a t e the accuracy of the s i m u l a t i o n of the f i e l d c u r r e n t , an attempt was made to d u p l i c a t e a f i e l d t e s t [54], In t h i s t e s t , the generator described i n example 2 was connected to a system as shown i n F i g . 23, and a three-phase f a u l t was a p p l i e d to the hig h s i d e of the step-up transformer. The network was simulated as coupled inductances w i t h the f o l l o w i n g parameters given i n [54]: zero sequence inductance L Q = 0.22 H, p o s i t i v e sequence inductance = 0.096 H. The i n i t i a l v o l t a g e at the generator terminals was 13.8/-30 0 kV (RMS, l i n e - t o - l i n e ) , and the i n i t i a l f i e l d c urrent was 620 A. The s w i t c h i n g sequence was as f o l l o w s [54]: Phase "b" at t = 0 s, phase "c" at t = 6 ms, and phase "a" at t = 20 ms. The simulated f i e l d current i s compared w i t h the measured f i e l d c urrent i n F i g . 24. (a - f i e l d t e s t , b — s i m u l a t i o n ) Examples 2 and 3 v e r i f y t h a t the generator model as w e l l as the numerical approach of method I give reasonably accurate r e s u l t s . I t i s , t h e r e f o r e , p o s s i b l e to proceed to a comparison of method I I w i t h method I . 66. 3 SINGLE PHASE UNITS 13.8/130.8 kV L =63.2mH INFINITE BUSBAR TRANSMISSION NETWORK L0 = 0.22H \ L1 = 0. 096 H FAULTED BUSBAR F i g . 23. System diagram. b- 10 a b TIME fs) 0 0.1 0.2 F i g . 24. Comparison of the simulated and measured f i e l d c u rrent. (5) Example 4 f o r T e s t i n g Method I I In t h i s example, a s i n g l e l i n e - t o - g r o u n d f a u l t was s t u d i e d . The system was simulated as shown i n F i g . 23. However, some changes were introduced. The i n f i n i t e busbar v o l t a g e was increased to 137.23/-20 kV (RMS, l i n e - t o - l i n e ) and the parameter X^" was changed from 0.29 p.u. to 0.18 p.u., which i s a more r e a l i s t i c value f o r a generator w i t h damper windings, than the previous one given i n [53]. The r e s u l t s were obtained i n three d i f f e r e n t ways: (a) S i m p l i f i e d generator model (voltage source E^^cos(tbt4p ) behind an impedance + J O J L ^ " ) ; (b) D e t a i l e d generator model i n t e r f a c e d w i t h method I (simultane-ous s o l u t i o n ) ; (c) D e t a i l e d generator model i n t e r f a c e d w i t h method I I (new tech-nique) . F i g . 25 compares the simulated current i n the f a u l t e d phase "a", and F i g . 26 compares the simulated f i e l d c u r r e n t . 25.0r -25.0\ 1 . 1 i i i 0.0 0.04 0-08 0.12 TIME (s) F i g . 25. Comparison of the simulated current i n the f a u l t e d phase "a". Figures 25 and 26 imply t h a t the s i m p l i f i e d g e n e r a t o r model i s r e a s o n a b l y a c c u r a t e f o r s h o r t time s t u d i e s , but no I n f o r m a t i o n can be o b t a i n e d on 6 8 . f l • • 1 , , , 0 0-04 0-08 0.12 TIME (s) F i g . 26. I d e n t i c a l r e s u l t s f o r f i e l d c u rrent. r o t o r c i r c u i t q u a n t i t i e s . The r e s u l t s obtained w i t h method I I are i n d i s t i n g u i s h a b l e from the r e s u l t s obtained w i t h model I . (6) Example 5 f o r Testing Method I I This example should provide a more severe t e s t f o r method I I , because i t a l s o includes t r a v e l l i n g wave e f f e c t s . The generator from example 4 was connected to a system, as shown i n F i g . 27. The transmission l i n e had the f o l l o w i n g parameters ("0" = zero sequence, "1" = p o s i t i v e sequence): 69. R = 0.2026 n/km , L =2.749 mH/km , C = 6.326 uF/km R = 0.0886 fi/km , L = 1.005 mH/km., C l = 11.408 yF/km and a length of 96.72 km. The length of the l i n e was chosen i n such a way as to match the p o s i t i v e sequence inductance of the transmission network from example 4. , 13.8/130.8 kV (3 single phase units ^-faulted busbar A<t . /TRANSMISSION LINE 1=4.83 JL-1 3.6° k A L=63.2mH infinite busbar-^ /seen from I side high) V=137.231^20° kV ' (line to line) F i g . 27. System diagram. Once more a single line-to-ground fault on the high side of the step-up transformer was studied. Fig. 28 compares the generator current in the faulted phase "a", calculated in three different ways: (a) Simplified generator model; (b) Detailed generator model interfaced with method I (simultane-ous solution); (c) Detailed generator model interfaced with method II (new technique). Again, the results from methods I and II are practically indistinguish-able . Fig. 28 may lead to the conclusion that the simple generator model i s as good as the detailed model. This i s only true, however, for 70. SIMPLIFIED MODEL 25.0, 12.5 •ic Uj Q: - 72.5 U -25.0 SIMULTANEOUS SOLUTION NEW TECHNIQUE 8.0 12.0 16.0 20.0 240 TIME (ms) Fig. 28. Comparison of the simulated current for faulted phase "a". some studies conducted over a very short time span where the flux decay does not play an important role. Fig. 29 compares the three-phase i n -stantaneous power at the generator terminals for the same case and shows that the simple model is much less adequate when power is measured. Finally, Fig. 30 shows the f i e l d current calculated with the two different interface techniques. As for the stator current, methods I and II give again results which are practically identical. The last two examples prove that the results obtained with method II agree to a high degree with results obtained with method I. This in turn proves the adequacy of method II. 71. SIMPLIFIED MODEL SIMULTANEOUS SOLUTION NEW TECHNIQUE TIME (ms) Fig. 29. Comparison of the simulated three-phase instantaneous power. 72. 600r u } - — 1 1 1 1 1 i i i 0 0.009 0.018 0-027 0.036 TIME (s) Fig. 30. Identical results for f i e l d current. The set of numerical examples would not be complete without a presentation of a case with a multi-mass mechanical system. A benchmark test case for such a system became available after completion of the thesis. It i s , therefore, not included in the main body of the text, but is added as Appendix 8. 73. 5. INITIAL CONDITIONS, DATA SCALING AND SATURATION 5.1 Calculation of the I n i t i a l Conditions for a Synchronous Generator The state of a synchronous generator is f u l l y described by the following variables: (1) a l l currents and voltages; (2) a l l angular displacements and speeds of the shaft system. The stator currents and voltages are normally obtained from a phasor solution of the entire system in which the generator i s represented as sinusoidal voltage or current sources. This solution usually considers the fundamental frequency only. The rotor circuit variables have to be found from the generator equations, which is straightforward for the case of balanced steady-state operation, but more complicated for the un-balanced case. In the latter case, harmonics exist not only in the rotor circuits but also in the stator circuits [8], [9], which leads to contradictions i f the total system was solved at fundamental frequency only. It is common practice, therefore, to assume a balanced steady-state operation for the generator. In this case, the damper currents are zero and the f i e l d current is constant. The currents and voltages w i l l vary sinusoidally in phase coordinates, but are constant in d,q,0-coordinates. It i s , therefore, easier to use d,q,0-coordinates for the calculation of the i n i t i a l conditions. The general voltage equations of the generator were defined in section 2.1 (eq. 21). They are shown again to aid understanding: v = -[L ]4— i - [R]i - [L']i (96) p dt —p —p p —p For the special case of balanced, steady-state conditions, where 1 = 0 , 74. they can be r e w r i t t e n as f o l l o w s : v = - [ R ] i - [ L * ] i ~p -p p -p (96a) With damper currents and zero sequence voltage and current being zero, (96a) can be reduced to the f o l l o w i n g s e t of three voltage equations: v, = -R i , - coL i d a d q q (97) and v = -R i + toL.i, + u)/-£ M i q a q d d v / 2 f f Vf. = "Vf (98) (99) The c o e f f i c i e n t i n f r o n t of M^ i n (98) r e s u l t s from the normalized form of Park's transformation m a t r i x [P] as shown i n (19). To f i n d the current i ^ and the r o t o r angular p o s i t i o n 3(0) or 6 ( 0 ) , i t i s necessary to r e l a t e (97-99) i n d,q,0-coordinates to the phasor diagram of the machine shown i n F i g . 31. DIRECT AXIS QUADRATURE AXIS R E F E R E N C E AXIS FOR NETWORK PHASOR SOLUTION F i g . 31. Phasor diagram of a synchronous generator f o r balanced, steady-state o peration (R and X X not to s c a l e ) . ct Q. CJ 75. From (97) and (98) i t follows that: V q + ^ d " - R a ( ± q + ^ " J X q ( ± q + ^ + ( X d " V *d + 6 q a ° 0 ) where \ - 5/l M f h For a balanced steady-state operation, d,q,0-coordinates are related to phasor values by: and i + j i d = v/J'I e j < S (102a) v q + j v d = /3 V e~3S (102b) where I and V are RMS positive sequence phasors. From (102a) and (102b) i t follows that (100) can be transformed to the re-ference frame for the network phasor solution, e e j 6 + (X , - X. ) i , e j 6 = V + R l + j X l (103) q d q d a J q where the phasors e e ^ and (X, - X )i, , e ^ l i e on the quadrature axis. q d q d Equation (103) allows, therefore, the determination of the angle 6(0), i f the phasors V and I are known. The l e f t hand side of (103) is not impor-tant in i t s value, but i t s position is that of the quadrature axis. The rest of the ele c t r i c a l variables can then easily be calcu-lated from the following relationships: i d = /3 |7|'sin(Y* - 6) (104a) i = /3 |I| . C O S ( Y ' - 6) (104b) and v d = /3 |v| sin(a' - 6) (104c) v q - /3 |v| cos(a' - 6) (104d) 76. F i n a l l y , the f i e l d current i ^ can be found from (98): v + R i - X,i, ± = _3 aq d d (105) l3 7 2 M f The i n i t i a l conditions of the e l e c t r i c part of the generator are then f u l l y defined. I t remains, therefore, to determine the i n i t i a l conditions for the mechanical variables of the generator. The mechanical equations of the generator were defined i n section 2.5 ((35)-(39)). However, the equation f o r mass i i s shown again to aid understanding: d2e. de , , J . - — \ + D . . — — + D . , .—(e. - e. .) + D. (e. - e.,.) + 1 dt 1 1 d t * l d t 1 1 _ ^ " " d t 1 1 + 1 K i - i > i < e i - 8 ± + i > + K i , i + i ( 9 i " e i + i > " T i ( 1 0 6 ) For steady-state conditions (106) can be s i m p l i f i e d as follows: de. D. . — - + K. .(6. - 6. ) + K. (6. - 9 ) = T. (106a) n d t I i - l ' i , i + l I i + l y I V ' dek The angular speed i s equal f o r a l l r o t a t i n g masses and can be found from the following r e l a t i o n s h i p : d 6 k 2 — — = u = u . £. for k = 1, ... N (107) dt m n ' where "subscript "m" denotes mechanical v a r i a b l e s , where n i s the number of poles i n the generator, and where CJ i s the angular frequency of the network. The i n i t i a l angular p o s i t i o n of the generator rotor can be cal c u l a t e d from the angle 5(0) as follows: 6 r = (6(0) + f)« f (108) where subscript " r " denotes generator rotor variables. Finally, the angle 6^ can be found from the angle 0^ : i-1 i-1 I T m3 " I D co m K i " l , i In a similar way, the angle can be found from: (109) N N I Tmj - I D a , K i , i + 1 (110) The sum of the applied mechanical torques Tmj must, of course, equal the sum of el e c t r i c a l and speed self-damping torques, so that there is zero accelerating torque i n i t i a l l y : N N N y T . = y Tej + y D..U> mi) ,L-. mi >, J . L , i i m J=l J j-1 J=l J J where subscript "e" denotes electromechanical torque. Calculation of the i n i t i a l angular displacement of the masses in the shaft system ends the process of i n i t i a l i z a t i o n of the generator variables, 5.2 Consistent Per Unit (p.u.) System and Conversion to Physical Units The choice of a consistent and simple p.u. system i s , in general, relatively easy. For rotating machinery, however, the situa-tion gets complicated i f more than one reference frame is used. A transformation from one reference frame to another may limit the freedom in the choice of base variables, i f symmetry of matrices is to be pre-served. This is precisely the case for a synchronous generator when unnormalized transformation matrices are used. 78. It i s common practice in the power industry to describe the generator in Park's d,q,0-coordinates, as explained in section 2.2. The conventional unnormalized transformation does not preserve the symmetry of the inductance matrix [L] of the generator [17]. The resulting asymmetrical matrix [L ] can be forced back to symmetry with a specific 3 p.u. system in which base power for rotor quantities is ^ - times base power for stator quantities [17], [55], A simpler approach, which does not require complicated scaling procedures to restore the symmetry, is presented in this chapter. As mentioned in section 2.2, the normalized transformation defined in (18) and (19) preserves the symmetry of the inductance matrix [L], i.e., the resulting matrix [L p] is symmetrical no matter which base values are chosen for stator and rotor circuits [17]. Then the conversion to p.u. values is a simple scaling problem with com-plete freedom in the choice of base values for yeach cir c u i t . Any linear electric network in steady-state operation can be described by one of the following nodal equations: or [Z . M . = V , (113) phys -^ phys -^ -phys where the subscript "phys" denotes physical values, and \ t Zphys ] = ^Phys 3" 1 ( 1 1 4 ) and where V , and I , are vectors of nodal voltages and currents —phys -phys injected into the nodes, in V and A respectively. The nodal impedance matrix [Z ^ ] is given in i t s own physical units, fi. In general, a system has more than one voltage level, with coupling through transformers. Therefore, different base voltages are normally 79. chosen for the node voltages, while base power i s normally the same f o r a l l nodes. I t can e a s i l y be shown that the following r e l a t i o n s h i p holds between the p h y s i c a l and the p.u. values [56]: [ Z P . U . ] " [ v b r l r W [ V " X ] <15> or where: [Z . ] = [V. ][Z ][V. ]• [S, ] 1 (116) phys b p.u. b b v 7 [V^] = diagonal matrix of base voltages, [S^] = diagonal matrix of base powers. Equations (115) and (116) are v a l i d f o r any set of base voltages and powers. However, symmetry of the matrix [Zp^yg] w i l l he preserved only i f there i s only one base power, which i s the normal p r a c t i c e i n power system analysis anyhow. I t i s customary to define the data of a synchronous generator i n a p.u. system based on i t s nameplate ratings. For general network studies, where each element has i t s own nameplate r a t i n g , a l l values must e i t h e r be converted to the same bases, or to p h y s i c a l q u a n t i t i e s . Per unit values o f f e r advantages i f a problem i s studied on a network analyzer, or on a d i g i t a l computer with fixed-point arithmetic, because i n both cases a l l values must be of a c e r t a i n order of magnitude.• This problem does not e x i s t i n computers with f l o a t i n g - p o i n t arithmetic, which i s the r u l e nowadays. Scaling (conversion from p.u. to p h y s i c a l values or v i c e versa) has no influence on the s o l u t i o n process, except f o r possible differences i n the accumulation of round-off e r r o r s . As a consequence, p r a c t i c a l l y i d e n t i c a l solutions (except f o r s l i g h t . d i f -ferences i n round-off errors) w i l l be obtained with p h y s i c a l quantities 80. and with p.u. quantities. The influence of scaling on round-off errors i s not easy to assess [57]. For a system of linear equations, however, the following statement can be proved [58]: "If scaling is done in such a way that only the exponent changes in floating-point numbers and i f the order of elimination is not changed, then the scaled (p.u.) equa-tions w i l l produce precisely the same significants in a l l answers and in a l l intermediate numbers". After careful examination of a l l advan-tages and disadvantages i t was decided to convert the machine data to physical units. Physical quantities are least confusing and assure con-sistency with the Transients Program. Conversion to physical units is done as follows: (a) Calculate a l l the required generator parameters in p.u. from the p.u. data based on nameplate ratings as provided by the manufacturer; (b) Multiply a l l elements of the matrices [L ] and [R] by p p.u. p.u. J the base impedance of the stator windings, which can be found from the following relationship: V 2 zsb = s?r <117> bb where: = base impedance of the stator; ^three-phase rated power of the stator in MVA, SSb "< i f the stator is onnected in wye. single-phase rated power of the stator in MVA, i f the stator is connected in delta. = line-to-line rated voltage of the stator in KV (RMS) for both wye and delta connections. 81. This operation transforms a l l the stator data to their o r i -ginal physical values. Rotor parameters, as found in step b w i l l be in physical values referred to the stator side of the generator. (c) Multiply a l l the parameters related to the rotor circuits 2 which l i e on the diagonal of the matrices [ Lpl a n d [R] by n » and those on the off-diagonal of [L^] by n, where: n = transformer ratio between the stator and the f i e l d . This number can be calculated in one of two ways: (1) from the physical value of the f i e l d resistance R^ , i f such i s known from measurements [22]: or: (2) from the open-circuit characteristic, as schematically shown in Fig. 32. P (118) A I R - G A P L I N E 1.0 Fig. 32. Typical open-circuit characteristic. From the generator equations derived in chapter 2.2, i t follows that the physical value of the field-to-stator mutual coupling is given by: n f J/f-Hj (119); where V = RMS value of line-to-line voltage found on the open-circuit Lib characteristic. From (119), the transformer ratio n can be defined as follows: M n = - (120) fp.u. Sb The procedure outlined above assumes that the original p.u. data was a l l based on the same base power, which is normally true for manufacturer's data. Then, the matrix [S^] is simply a unit matrix premultiplied by a scalar Sg^ defined in (117) . This procedure also assumes that there are only two base voltages, one for the stator, and the other for the f i e l d and damper windings. The latter assumption can be j u s t i f i e d as follows: Since the damper windings are hypothetical windings, for an inter-connected arrangement of many damper bars, any transformer ratio to them can be assumed. It i s , therefore, possible to use the same ratio as that from the stator to the f i e l d without loss of generality. This specific conversion procedure is very simple to use and requires only the available standard data. A similar approach based upon d i f -ferent reasoning has been suggested in [59]. 83. 5.3N S a t u r a t i o n i n the Steady-State Operation of a Synchronous Generator S a t u r a t i o n may have an impact (sometimes a s i g n i f i c a n t one) on power system t r a n s i e n t s t a b i l i t y and steady-state s t a b i l i t y c a l c u l a -t i o n s , as w e l l as on r e a l and r e a c t i v e power flow c a l c u l a t i o n s [29], [60]. As i s w e l l known i n p r a c t i c e , i t also i n f l u e n c e s the c a l c u l a t i o n s of electromagnetic t r a n s i e n t s i n power systems [61]. On the other hand, a n o n l i n e a r r e l a t i o n between f l u x and current w i l l , i f t r e a t e d r i g o r -o u s l y , complicate the s o l u t i o n process s i g n i f i c a n t l y . The s o l u t i o n of l a r g e n o n l i n e a r systems becomes then very expensive and time consuming. An approximate treatment of s a t u r a t i o n e f f e c t s i s , t h e r e f o r e , commonly accepted. The treatment of s a t u r a t i o n e f f e c t s i n steady-state operation d i f f e r s from the approach needed f o r t r a n s i e n t s i m u l a t i o n s . Before t r e a t i n g the l a t t e r case i n s e c t i o n 5.4, i t seems appropriate to present a short review of the e x i s t i n g approaches used i n s t a b i l i t y s t u d i e s , which f a l l i n t o the category of steady-state phasor s o l u t i o n s . One of the e a r l i e s t approaches towards s a t u r a t i o n i s t o be found i n [18], [62], where on an e m p i r i c a l b a s i s , i t was suggested to c a l c u l a t e the values of the s a t u r a t e d reactances by m u l t i p l y i n g the unsaturated values by a constant F = 0.88. This simple approach i s c l e a r l y not accurate enough, s i n c e the s a t u r a t i o n e f f e c t s vary w i t h the type of generator and i t s l o a d i n g c o n d i t i o n s [63]. More recent' approaches can be d i v i d e d i n t o the f o l l o w i n g b a s i c groups: (1) The degree of s a t u r a t i o n i s a f u n c t i o n of the t o t a l f l u x f~2 2 ty = / ij;^ + ij; . There i s , t h e r e f o r e , only one s a t u r a t i o n c o e f f i c i e n t f o r the t o t a l f l u x [64-66]; (2) The degree of s a t u r a t i o n i n each a x i s i s p r o p o r t i o n a l to the components of the v o l t a g e source behind P o t i e r (or leakage) 84. reactance. There are, t h e r e f o r e , two separate s a t u r a t i o n c o e f f i c i e n t s , one f o r each a x i s [24], [67], A number of authors use f l u x p l o t s to determine the s a t u r a -t i o n e f f e c t s [60], [68], The r e s u l t s published i n [68] seem to favour approach 1. I t i s also worth mentioning that s a t u r a t i o n data f o r a medium-sized generator published i n [69] are clos e to the value of the e m p i r i c a l c o e f f i c i e n t F suggested i n [18]. Other authors, however, suggest d i f f e r e n t s a t u r a t i o n c h a r a c t e r i s t i c s f o r the d- and q-axes. To sum i t up, i t i s not yet known which procedure i s more accurate. The main problem l i e s i n the u n a v a i l a b i l i t y of accurate data [29]. A p o s s i -b l e s o l u t i o n of t h i s problem could come from measuring the s a t u r a t i o n e f f e c t s d i r e c t l y i n phase coordinates [12], [22], but more research i s needed [60]. The n o n l i n e a r f l u x - c u r r e n t c h a r a c t e r i s t i c caused by s a t u r a -t i o n i m p l i e s t h a t i t i s no longer p o s s i b l e to use, i n a s t r a i g h t f o r w a r d way, phasor s o l u t i o n s i n the c a l c u l a t i o n of the steady-state c o n d i t i o n s . To get around t h i s problem, an "eq u i v a l e n t l i n e a r machine", which i s e x a c t l y v a l i d i n one p a r t i c u l a r operating p o i n t and approximately v a l i d i n the neighbourhood of that p o i n t , i s introduced [70]. The o b j e c t i v e i s to l i n e a r i z e the problem by parameter m o d i f i c a t i o n . This i s achieved by r e p l a c i n g the n o n l i n e a r c h a r a c t e r i s t i c by a l i n e a r curve through the ope r a t i n g p o i n t and the o r i g i n , as shown i n F i g . 33. This approxima-t i o n i s , of course, v a l i d only i n the immediate v i c i n i t y of the opera-t i n g p o i n t . The f o l l o w i n g procedure was developed f o r the i n c l u s i o n of s a t u r a t i o n e f f e c t s i n t o the c a l c u l a t i o n of the i n i t i a l c o n d i t i o n s of a generator: 10 APPROXIMATE / L I N E A R ^ J — - C U R V E OPERATING POINT i f F i g . 33. L i n e a r i z a t i o n through the o r i g i n . (1) C a l c u l a t e the i n i t i a l c o n d i t i o n s assuming that the generator i s unsaturated, as e x p l a i n e d i n s e c t i o n 5.1. (2) I f the generator operates i n the s a t u r a t e d r e g i o n , the c a l c u -l a t i o n s are repeated w i t h the unsaturated parameters replaced by t h e i r s a t u r a t e d equivalent values. These values can be found from the slope of the approximate s t r a i g h t l i n e shown i n F i g . 33. The process i s repeated u n t i l i t converges. The steady-state conditions of the network found w i t h phasor s o l u t i o n techniques do not contain harmonics. By s i m u l a t i n g the problem i n steady-state (no f a u l t applied) as a t r a n s i e n t case f o r a p e r i o d of a few c y c l e s , s t a r t i n g from the i n i t i a l c o n d i t i o n s w i t h the l i n e a r i z e d generator, a new steady-state w i t h harmonics should be reached i f s a t u r a t i o n i n generators and transformers i s modelled i n the Transients Program [45]. This approach worked q u i t e w e l l i n a case where the transformer s a t u r a t i o n generated harmonics [71], but i t has not yet been t e s t e d f o r generator s a t u r a t i o n . A s i m i l a r approach could be used f o r 86. unbalanced cases, s i n c e the i n i t i a l conditions are c a l c u l a t e d f o r p o s i -t i v e sequence currents o n l y , as explained i n s e c t i o n 5.1. 5.4 D e f i n i t i o n s of S a t u r a t i o n i n the Simulation of Electromagnetic Transients The a n a l y s i s of t r a n s i e n t performance of synchronous generators w i t h constant inductances may l e a d to s e r i o u s e r r o r s both i n form and magnitude of currents and voltages [61], Even the use of the term "inductance" may be m i s l e a d i n g , s i n c e i t i s based upon the assumption of l i n e a r i t y , i . e . , i t i s no longer true that ty = L*MMF. I t should, r a t h e r , be s a i d that ty i s a n o n l i n e a r f u n c t i o n of MMF. For example, [22] i n t r o -duces two types of inductances: (a) secant inductance, defined by t o t a l f l u x per u n i t c u r r e n t . (b) incremental inductance, defined by the r a t e of the change of f l u x l i n k a g e w i t h respect to current. I t i s proposed to base the a n a l y s i s of s a t u r a t i o n e f f e c t s upon the f o l l o w i n g assumptions: (1) In any reference frame, the generator f l u x e s can be represented as f o l l o w s : 1 = i L + (121) where: 3^= v e c t o r of f l u x e s r e l a t e d to leakage inductances, uneffected by s a t u r a t i o n ; v e c t o r of f l u x e s r e l a t e d to mutual inductances, subjects to s a t u r a t i o n e f f e c t s . Only the l a t t e r f l u x e s w i l l be considered i n the f o l l o w i n g a n a l y s i s . (2) The degree of s a t u r a t i o n i s a f u n c t i o n of the MMF, which i n turn is a function of the total unsaturated flux ty calculated u along the airgap line. (3) The saturation effects are equal on both.axes, i.e., there is only one saturation coefficient. (4) The distortion of any airgap flux waves does not effect the unsaturated inductance values or destroy the sinusoidal varia-tions assumed for rotor and stator inductances. (5) Hysteresis and eddy current losses are neglected, as is usual in power transformer modelling where i t has normally l i t t l e influence on the results [72], Assumption 2 implies the knowledge of the dependence between the instan-taneous flux and the excitation current. The only available data, how-ever, consists of the open circuit characteristic (terminal voltage as a function of the excitation current) shown schematically in Fig. 32. The converted curve (flux versus current) has the same form as the o r i -ginal curve. A short proof is given in Appendix [7], e.g., straight-line segments, exponential or quadratic curves. It was decided to adopt a two straight-lines approximation, due to i t s simpli-city, but the actual number of segments does not change the method of analysis. It can be increased, i f so required by the shape of the flux-current curve. The total (mutual) flux of the machine may, there-fore, be described as follows (subscript "u" denotes an unsaturated value, subscript "m" dropped to simplify notation): The resulting curve can be approximated in a number of ways, unsaturated region (122) m ty + a T i t saturated region 88. where: a and m = constants resulting from the two straight-lines approximation of the saturation curve, or from an approximation with more than two linear segments. The total unsaturated airgap flux, on the other hand, may be described by the following equation: ' ty - /ty] + ty2 (123) u du qu where subscripts "d" and "q" denote the direct and quadrature axis values, respectively. Equations (122) and (123) imply that there is one saturation effect for the total flux, rather than two separate effects, one for each axis. This situation is show schematically in Fig. 34. du UJ Uj rds rdu Fig. 34. Schematic representation of saturation effects. From Fig. 34 and equations (122) and (123) follows that the fluxes in the saturated region can be described as follows (subscript "s" denotes saturated values): i i , = n i ' i l i , + a cos 3 (124) rds du ij; = m«ij; + a sin qs qu where: and cos 3 = Kdu sin ty = -3S. (125) (126) (127) As already mentioned, the constants m and a result from the straight-line approximation of the flux-current characteristic. This situation is illustrated schematically in Fig. 35. Fig. 35. Straight-line approximation of the flux-current characteristic. From basic analytic geometry follows that m = m, (128) The constant a can be found from the conditions in the "knee-point" ( i j ; , i £ ) . In this point both fluxes (saturated and unsaturated) must be equal, i.e., m i = m i + a 1 c 2 c (129) 9 0 . From ( 1 2 9 ) follows that: a = (m: - m 2 ) i c ( 1 3 0 ) Equations ( 1 2 4 ) - ( 1 3 0 ) describe f u l l y the saturated fluxes in a synchron-ous generator. They can be easily modified to accommodate additional data, i f such i s available. It i s , for example, possible to create two sets of constants m and a, one for each axis. This would allow each axis to have i t s own saturation coefficient. 5.5 Implementation in the Transients Program The generator equations in d,q,0-coordinates can be rewritten into the following form (details were given in section 2 . 2 ) : ( 1 3 1 ) where the matrix [A] is defined as follows: 0 CO 0 1 -CO 0 o 1 1 0 ( 1 3 2 ) 0 0 1 0 _ _|_ _ 0 1 0 [A] = [P]~[P] 1 Only the two last terms of ( 1 3 1 ) are subject to saturation influence. Their different physical nature results in two different implementation procedures, one for each of them. Both procedures, however, are based upon the following common assumption: The generator does not change i t s saturation status during one time-step, i.e., i f the generator was saturated at the beginning of a time-step, i t remains saturated at i t s end. The saturable transformer voltages are described by the following equation: 91, v = " 4r JL (133) —pm dt -^-pm Equation (133) i m p l i e s that p h y s i c a l l y only the incremental f l u x e s are of importance. Since the generators equations are solved w i t h the t r a p e z o i d a l r u l e of i n t e g r a t i o n , (133) i s then transformed to the f o l l o w i n g form ( s u b s c r i p t "p" dropped to s i m p l i f y n o t a t i o n ) : i L ( t ) = j ; ( t - A t ) - ~ ( v ( t ) + v ( t - A t ) ) (134) " T i l 2. m m S u b s t i t u t i o n of (122) i n t o (134) y i e l d s the f o l l o w i n g expression f o r the i t h component of the v e c t o r ^ ( t ) ( s u b s c r i p t "m" dropped to s i m p l i f y n o t a t i o n ) : m ty. (t) + a = m ty. ( t - A t ) + a - - ^ ( v . ( t ) + v. ( t - A t ) ) (135) I U i u 2 l l Simple rearrangements y i e l d the f o l l o w i n g r e s u l t : v i ( t ) = ^ i r ^ i u ^ _ * i u ( t - A t ) ) - v i ( t _ A t ) < 1 3 6> Equation (136) provides the means f o r i n c l u d i n g the s a t u r a t i o n e f f e c t s i n the transformer voltages of a generator. I t i s simply enough to m u l t i p l y a l l the mutual inductances by the constant m. I f the parameters of r o t o r c i r c u i t s of the generator are not r e f e r r e d to the s t a t o r s i d e , i t i s necessary to introduce the transformer r a t i o i n t o the constant m. The sa t u r a b l e terms r e l a t e d to the speed voltages appear only i n the f o l l o w i n g two equations: and v, = -co ty (137) dm rqm v ' v = +co ty. (138) qm Tdm v ' where s u b s c r i p t s "d" and "q" denote the d i r e c t and quadrature a x i s , r e s p e c t i v e l y , and s u b s c r i p t "m" denotes terms r e l a t e d to mutual induc-tances. Equations (137) and (138) imply that i t i s necessary to con-s i d e r the e n t i r e f l u x e s ty, and ty , since (137) and (138) represent dm qm a l g e b r a i c r e l a t i o n s h i p s , r a t h e r than d i f f e r e n t i a l r e l a t i o n s h i p s . I f a l l r o t o r c i r c u i t s are converted to the s t a t o r s i d e , the f o l l o w i n g r e l a t i o n s h i p s can be obtained f o r the saturate d f l u x e s ( s u b s c r i p t "m" dropped to s i m p l i f y n o t a t i o n ) : ty, = m M , ( i , + i £ + i ) + a cos 6 (139) ds r d r D and ty = m M ( i + i + i . ) + a s i n 0 (140) qs q q g Q I t i s , t h e r e f o r e , enough to s u b s t i t u t e (139) and (140) f o r the f l u x e s r e l a t e d to the sa t u r a b l e inductances. The procedure o u t l i n e d above i s very f l e x i b l e and can e a s i l y accommodate a d d i t i o n a l data, as i t becomes a v a i l a b l e . I t i s not t i e d e n t i r e l y t o the theory o u t l i n e d i n the previous chapter, and i t can be adapted to accept any a v a i l a b l e r e p r e s e n t a t i o n of s a t u r a t i o n e f f e c t s . 93. 6. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH A range of problems related to interfacing generator models with an Electromagnetic Transients Program has been described. The main topics covered in this thesis were: (a) creation of an adequate generator model (b) development of two alternative interfacing techniques (c) modelling of saturation effects in the generator. The validity of the generator model has been verified in a number of test studies, which also included comparisons with f i e l d test results. Good agreement was achieved, and the results obtained with the two interfacing methods gave practically identical answers which proves their v a l i d i t y . The discussion of saturation effects presents only a f i r s t attempt i n this area. The suggested solution methods should be tested in practice to establish their validity. Similarly, the proposed method II of interfacing generator models with the Transients Program should undergo further tests before completely replacing method I with i t . Further possible areas of additional investigations could include such topics as: (a) inclusion of space harmonics in the f i e l d distribution (b) improvements in the calculation of i n i t i a l conditions to allow the i n i t i a l i z a t i o n from unbalanced conditions. This work should, therefore, be understood as a f i r s t step, only, in the modelling of synchronous generators for electromagnetic transient phenomena. 94. REFERENCES [1] S. Meyer and H.W. Dommel, "Numerical M o d e l l i n g of Frequency-Dependent Transmission-Line Parameters i n an Electromagnetic Trans-i e n t s Program", IEEE Trans. Power App. Syst., v o l . PAS-93, Sept./ Oct. 1974, pp. 1401-1409. [2] V. Brandwajn and H.W. Dommel, " I n t e r f a c i n g Generator Models w i t h an Electromagnetic Transients Program", Paper A 76359-0, presented at the IEEE PES Summer Meeting, P o r t l a n d , Ore., J u l y 19-23, 1976. [3] K. Carlson, E.H. L e n f e s t , and J . J . La F o r e s t , "MANTRAP, Machine and Network Transients Program", Paper V-D, presented at the 9th PICA Conference, New Orleans, La., June 2-4, 1975. [4] M.C. H a l l , "Machine Model f o r Transients Program", I n t e r n a l Memor-andum, Southern C a l i f o r n i a Edison Company, Rosemead, Ca., 1975. [5] IEEE Subsynchronous Resonance Task Force Report, " F i r s t Benchmark Model f o r Computer Simulation of Subsynchronous Resonance", Paper F 77102-7, presented at the IEEE PES Winter Meeting, New York, N.Y., Jan. 30-Feb. 4, 1977. [6] P. Eykhoff, "System I d e n t i f i c a t i o n , Parameter and State E s t i m a t i o n " , W i l e y - I n t e r s c i e n c e , N.Y., 1974. [7] R.H. Park, "Two Reaction Theory of Synchronous Machines", AIEE Trans., v o l . 48, p t . I , J u l y 1929, pp. 716-730. [8] G. M u l l e r , "Theorie der rotierenden Maschine", L e i p z i g , 1968. [9] B. Adkins and R.G. Harley, "The General Theory of A l t e r n a t i n g Cur-rent Machines: A p p l i c a t i o n to P r a c t i c a l Problems", Chapman and H a l l , London, 1975. [10] IEEE Committee Report, "Recommended Phasor Diagram f o r Synchronous Machines", IEEE Trans. Power App. Syst., v o l . PAS-88, Nov. 1969, pp. 1593-1610. [11] B. Adkins, "Tr a n s i e n t Theory of Synchronous Generators Connected to Power Systems", Proc. IEE, v o l . 98, No. 2, 1951, pp. 510-523. [12] M. Ra f i a n and M.A. Laughton, "Determination of Synchronous-Machine Phase-Co-Ordinate Parameters", Proc. IEE, v o l . 123(8), Aug. 1976, pp. 818-824. [13], Y. Takeda and B. Adkins, "Determination of Synchronous-Machine Parameters A l l o w i n g f o r Unequal Mutual Inductances", Proc. IEE, v o l . 121(12), Dec. 1974, pp. 1501-1504. [.14] I.M. Canay, "Causes of Discrepancies on C a l c u l a t i o n of Rotor Quan-t i t i e s and Exact E q u i v a l e n t Diagrams of the Synchronous Machine", IEEE Trans. Power App. Syst., v o l . PAS-88, J u l y 1969, pp. 1114-1120. 95. [15] H.Spath, "Der E i n f l u s s von Kopplungen zwischen dem Nullsystem und der d- order q-Achse auf das B e t r i e b von Synchron-Schenkelpolmas-chinen", ETZ-A, Bd 93, 1972, pp. 712-714. ( [16] E.W. Kimbark, "Power System S t a b i l i t y : Synchronous Machines", New York: Dover, 1968 ( r e p r i n t of 1956 ed.). [17] M.R. H a r r i s , P.J. Lawrenson, and J.M. Stevenson, "Per-Unit Systems w i t h S p e c i a l Reference to E l e c t r i c a l Machines", Cambridge U n i v e r s i t y P r ess, Cambridge, 1970. [18] L.A. K i l g o r e , " C a l c u l a t i o n of Synchronous Machine Constants", Trans. AIEE, v o l . 50, 1931, pp. 1201-1214. [19] C. Concordia, "Synchronous Machines: Theory and Performance", John Wiley and Sons, New York, 1951. [20] Y. Yu and H.A.M. Moussa, "Experimental Determination of Exact E q u i -v a l e n t C i r c u i t Parameters of Synchronous Machines", IEEE Trans. Power App. S y s t . , v o l . PAS-90, Nov./Dec. 1971, pp. 2555-2560. [21] G. Shackshaft, "New Approach f o r the Determination of Synchronous-Machine Parameters from Tests", Proc. IEE, v o l . 121(11), Nov. 1974, pp. 1385-1392. [22] L.A. Snider and I.R. Smith, "Measurement of Inductance C o e f f i c i e n t s of Saturated Synchronous Machine", Proc. IEE, v o l . 119(5), May 1972, pp. 597-601. [23] G. Manchur e t . a l . , "Generator Models E s t a b l i s h e d by Frequency Response Tests on a 555 MVA Machine", S t a b i l i t y of Large E l e c t r i c Power Systems, IEEE P r e s s , New York, 1974, pp. 37-44. [24] P. Dandeno, R.L. Hauth, and R.P. Schulz, " E f f e c t s of Synchronous Machine M o d e l l i n g i n Large Systems S t u d i e s " , i b i d , pp. 28-36. [25] V. Gourishankar, " E l e c t r o m a g n e t i c a l Energy Conversion", I n t e r n a -t i o n a l Textbook Company, Seranton, 1975. [26] L.A. Soderquist, "Fundamentals of System S t a b i l i t y " , The Role of Primes Movers i n Systems S t a b i l i t y , IEEE T u t o r i a l Course, New York, N.Y. , 1970, pp. 7-15. [27] R.T. B y e r l y , "Power Systems S t a b i l i t y - E f f e c t s of C o n t r o l System Performance", i b i d , pp. 57-66. [28] P. Subramian and P. M a l i k , " D i g i t a l S i m u l a t i o n of a Synchronous Generator i n Direct-Phase Quantities", Proc. IEE, v o l . 118(1), Jan. 1971, pp. 153-160. [29] C. Concordia, D i s c u s s i o n to G. Shackshaft and R. N e i l s o n , "Results of S t a b i l i t y Tests on an Underexcited 120 MW Generator", Proc. IEE, v o l . 119, 1972, pp. 1487-1494. 96. [30] F.H. B r a n i n , J r . , "Computer Methods of Network A n a l y s i s " , Proc. IEEE, v o l . 55, Nov. 1967, pp. 1787-1801. [31] C.W. Gear, "Numerical I n i t i a l Value Problems i n Ordinary D i f f e r e n -t i a l Equations", P r e n t i c e - H a l l , Englewood C l i f f s , N.J., 1971. [32] G. Gross and A.R. Bergen, "A Class of New M u l t i s t e p I n t e g r a t i o n Algorithms f o r the Computation of Power System Dynamical Response", IEEE Trans. Power App. Syst., v o l . PAS-91, Jan./Feb. 1977, pp. 293-306. [33] C.W. Gear, "The Automatic I n t e g r a t i o n of Large Systems of Ordinary D i f f e r e n t i a l Equations", The Digest Record of the J o i n t Conference on Mathematical and Computer Aids to Design, Anaheim, Ca., Oct. 27-31, 1969, pp. 27-58. [34] H.W. Dommel and W.S. Meyer, "Computation of Electromagnetic Tran-s i e n t s " , Proc. IEEE, v o l . 62, J u l y 1974, pp. 983-993. [35] G.G. D a h l q u i s t , "A S p e c i a l S t a b i l i t y Problem f o r L i n e a r M u l t i s t e p Methods", Nordisk T i d s k r i f t f o r Informations-Behandling, BIT ( 3 ) , 1963, pp. 27-43. [36] S.N. Talukdar, " I t e r a t i v e M u l t i s t e p Methods f o r Transient S t a b i l i t y S t u d i e s " , IEEE Trans. Power App. Syst., v o l . PAS-90, Jan./Feb. 1971, pp. 96-102. r-^ji A. Semlyen and A. Dabuleanu, "A System Approach to Accurate Switch-i n g Transient C a l c u l a t i o n s Based on State V a r i a b l e Component Model-l i n g " , IEEE Trans. Power App. Syst., v o l . PAS-94, M a r c h / A p r i l 1975, pp. 572-578. [38] S.C. T r i p a t h y and N.D. Rao, "A-Stable Numerical I n t e g r a t i o n Method f o r Transmission System T r a n s i e n t s " , Paper F 77 062-3 presented at the IEEE PES Winter Meeting, New York, N.Y., Jan. 30-Feb. 4, 1977. [39] E.J. Davison, "A High-order Crank-Nicholson Technique f o r S o l v i n g D i f f e r e n t i a l Equations", The Computer J o u r n a l , v o l . 10(2), Aug. 1967, pp. 195-197. [40] H.W. Dommel, " D i g i t a l Computer S o l u t i o n s of Electromagnetic Tran-s i e n t s i n S i n g l e - and Multiphase Networks", IEEE Trans. Power App. Sys t . , v o l . PAS-88, Apr. 1969, pp. 388-399. [41] 0.1. E l g e r d , " E l e c t r i c Energy Systems Theory: An I n t r o d u c t i o n " , McGraw-Hill, New York, N.Y., 1971. [42] R.W. Hamming, "Numerical Methods f o r S c i e n t i s t s and Engineers", I n t e r n a t i o n a l Student E d i t i o n , McGraw-Hill, Tokyo, 1962. [43] C.W. Gear, "The Automatic I n t e g r a t i o n of S t i f f Ordinary D i f f e r e n -t i a l Equations", Proceeding IFIPS Conference, Edinburgh, 1968, pp. 187-193. 97. [44] W.S. Meyer, Transients Program Memorandum, v o l . UJV, March 1974. [45] H.W. Dommel, "Nonlinear and Time-Varying Elements i n D i g i t a l Simu-l a t i o n of Electromagnetic T r a n s i e n t s " , IEEE Trans. Power App. Syst., v o l . PAS-90, Nov./Dec. 1971, pp. 2561-2567. [46] H.J. C a r l i n and A.B. Giordano, "Network Theory: An I n t r o d u c t i o n to R e c i p r o c a l and N o n r e c i p r o c a l C i r c u i t s " , P r e n t i c e - H a l l , Engle-wood C l i f f s , N.J., 1964. [47] D.A. Calahan, "Numerical Considerations f o r Implementation of Non-l i n e a r Transient C i r c u i t A n a l y s i s Program", IEEE Trans, on C i r c u i t Theory, v o l . CT-18(1), Jan. 1971, pp. 66-73. [48] H.W. Dommel, "A Method f o r S o l v i n g Transient Phenomena i n Multiphase Systems", Report 5.8, presented at 2nd Power System Computation Conference, Stockholm, 1966. See a l s o E l e c t r i c Research C o u n c i l , Transmission Line Reference Book: 345 KV and Above, E l e c t r i c Power Research I n s t i t u t e , Palo A l t o , Ca., 1975, pp. 292-293. [49] L. Dube and H.W. Dommel, "Si m u l a t i o n of C o n t r o l Systems i n an E l e c -tromagnetic Transients Program w i t h TACS", Paper submitted f o r PICA Conference, Toronto, Ont., May 24-27, 1977. [50] E. Bodewig, " M a t r i x C a l c u l u s " , 2nd e d i t i o n , North-Holland, 1959. [51] S.N. Talukdar, "METAP - A Modular and Expandable Program f o r Simu-l a t i n g Power System T r a n s i e n t s " , IEEE Trans. Power App. Syst., v o l . PAS-95(6), Nov./Dec. 1976, pp. 1882-1891. [52] H.W. Dommel and N. Sato, "Fast T r a n s i e n t S t a b i l i t y S o l u t i o n s " , IEEE Trans. Power App. Syst., v o l . PAS-91, July/Aug. 1972, pp. 1643-1650. [53] A. Charlton and G. Shackshaft, "Comparison of Accuracy of Methods f o r Studying S t a b i l i t y , N o r t h f l e e t E x e r c i s e " , E l e c t r a , No. 23, J u l y 1972, pp. 9-49. [54] G. Shackshaft and R. N e i l s o n , "Results of S t a b i l i t y Tests on an Under-excited 120 MW Generator", C e n t r a l E l e c t r i c i t y Generating Board, J u l y 1971. [55] A.W. Rankin, "Per-Unit Impedances of Synchronous Machines", AIEE T r a n s a c t i o n s , v o l . 64, Aug. 1945, pp. 569-572 and pp. 839-841. [56] H.W. Dommel, "Notes on Power System A n a l y s i s " , U n i v e r s i t y of B r i t i s h Columbia, Vancouver, B.C., 1974. [57] S.H. C r a n d a l l , "Engineering A n a l y s i s - A Survey of Numerical Proce-dures", McGraw-Hill, New York, N.Y., 1956, p. 37. [58] G.E. Forsythe and C B . Moler, "Computer S o l u t i o n of L i n e a r Alge-b r a i c Systems", P r e n t i c e - H a l l , Englewood C l i f f s , N.J., 1967, p. 39. 98. D.C. White and H.H. Woodson, "Electromechanical Energy Conversion", J . Wiley, New York, N.Y., 1959, p. 518. N.H. Demerdash .and H.B. Hamilton, "A S i m p l i f i e d Approach t o Deter-mination of Saturated Synchronous Reactances of Large Turbogenera-t o r s under Load", IEEE Trans. Power App. Syst., v o l . PAS-95(2), M a r c h / A p r i l 1976, pp. 560-569. R. Rudenberg, "Transient Performance of E l e c t r i c Power Systems, Phenomena i n Lumped Networks", MIT Pr e s s , Cambridge, Mass., 1970. S.H. Wright, "Determination of Synchronous Machine Constants by Test", Trans. AIEE, v o l . 50, Dec. 1931, pp. 1331. B. Adkins, D i s c u s s i o n to [60]. G. Shackshaft, "General-Purpose Turbo-Alternator Model", Proc. IEE, v o l . 110(4), A p r i l 1963, pp. 703-713. W.D. Humpage and T.N. Saha, "Digital-Computer Methods i n Dynamic-Response Analyses of Turbogenerator U n i t s " , Proc. IEE, v o l . 114(8), Aug. 1967, pp. 1115-1130. J.D. Schulz, W.D. Jones and D.N. Ewart, "Dynamic Models of Turbine Generators Derived from S o l i d Rotor Equivalent C i r c u i t s " , IEEE Trans. Power App. Syst., v o l . 92(3), May/June 1973, pp. 926-933. D. W. O l i v e , "New Techniques f o r the C a l c u l a t i o n of Dynamic S t a b i -l i t y " , IEEE Trans. Power App. Sy s t . , v o l . PAS-85, J u l y 1966, pp. 767-777. E. F. Fuchs and E.A. E r d e l y i , "Determination of Waterwheel A l t e r n a t o r Steady State Reactances from F l u x P l o t s " , IEEE Trans. Power App. Syst., v o l . PAS-91, 1972, pp. 1795-1802. E.F. Fuchs, D i s c u s s i o n to [60]. B. Adkins, D i s c u s s i o n to [60]. H. W. Dommel and W.S. Meyer,"Computation of Electromagnetic Tran-s i e n t s " , Proc. IEEE, v o l . 62, J u l y 1974, pp. 983-993. Cigre Working Group 13.05, "Modelling of s a t u r a b l e Elements f o r Transient Overvoltages S t u d i e s " , to be pu b l i s h e d i n E l e c t r a . "Test Procedures f o r Synchronous Machines", IEEE P u b l i c a t i o n No. 115, 1965. B.C. Kuo, "Discrete-Data C o n t r o l Systems", P r e n t i c e - H a l l , Englewood C l i f f s , N.J., 1970, pp. 336-344. B. Noble, "Applied L i n e a r Algebra", P r e n t i c e - H a l l , Englewood C l i f f , N.J., 1969. 99. APPENDIX 1 DEFINITIONS OF L' L" T' , T" AND L', L", T» , AND T" d' d' do' do q' q' qo' go Only the definitions of direct axis quantities w i l l be derived. The quadrature axis quantities can be obtained from the definitions of the direct axis by replacing the subscripts "d", " f " , and "D" with "q", "g", and "Q", and by replacing the voltage v^ with 0, i.e., v^ = 0, since the g-winding is permanently short-circuited. The quantities L^, 1/j, T ^ q and are equivalent parameters which are only defined for transient conditions following a disturbance. According to the IEEE (ANSI) standards [73] in the U.S.A. and similar standards elsewhere, a simultaneous three-phase short-circuit i s used to measure these quantities. It can be simulated with (21)-(23) by apply-ing voltages to the generator terminals which are equal in magnitude and of opposite sign to those existing in the balanced steady-state atrrated speed prior to the fault. 1. Subtransient Inductance and Open-Circuit Time Constants Immediately after the disturbance ( f i r s t few cycles), there w i l l be currents flowing i n the damper windings. The following assump-tions can be made during that short period: (a) no voltage regulator action yet, i.e., v^ = const (1-1) (b) constant flux linkages in the rotor circuits, i.e., \pf = const (1-2) ijjjj = const (1-3) The last assumption i s equivalent to neglecting the resistances in the f i e l d and damper windings, which cause a slow decay in 100. the fluxes. Immediately after the disturbance, this decay is negligible. Equations(22) can be used to express the currents i ^ and i ^ as functions of these constant flux linkages: -M, D L f L D " M f * f - f ¥ d *D-2Vd (1-4) Substitution of the above expressions into the equation for the flux in the d-axis yields: (1-5) 3 2 L f + LD " 2 M f ^ = (I, " 4 - i 2 T ± ) i d + K ^ f + K ^ D "d ^"d 2 " f M2 f D - M f where K l = / f M f ( L D - Mf) L f V 4 (1-6) and K =' 2 2 Mf(Lf - V L fV 4 (1-7) The sub transient inductance L" relates stator flux \b, to stator currents d d i ^ . It i s , therefore, defined as follows: T II L 3 „2 f D 2 M f L d " L d - 2 M f — -T f D " " f d-8) The open-circuit time constants define the decay of the rotor fluxes after the disturbance. They can be found from the two differential equations for rotor fluxes. From (21) follows that dt % _dt J R, 0 • + "V f i 0 D (1-9) 101. Substitution of (1-4) into (1^ -9) with i d = 0 (open-circuit) yields: di|>f d t -dt L f L D - 4 R f L D " R f M f *f - v f • + 0 (1-10) The open-circuit time constants are the reciprocals of the eigenvalues of the matrix in (1-10) and they can be found from the following equation: where a T M a T J a a = L f L D - Mf Solving (1-11) for T yields the following result: Tdo rptt do 1 , L f . \ , 1 I S V 2 . , « J 2 / VR, (1-11) (1-12) (1-13) with the positive sign of the root for T d o » and the negative sign for r p l l do' 2. Transient Inductance L' d After elapse of a few cycles, i t can be assumed that the damper winding current has already died out, i.e., 1 D = 0 (1-14) Therefore, only the f i e l d winding with an unchanged flux \> has to be considered. From (22) i t follows that under these conditions the f i e l d current i s given as: 2 " f f r f L, (1-15) ) 102. Substitution of (1-15) into the flux equation for the direct axis yields: 2 M p— *d= < L d " f ^ *<! V* Mf *f ( 1 _ 1 6 ) The transient inductance i s , therefore, defined as follows: M2 L d = L d " I if t1"17) The quadrature axis quantities are defined in an analogous way with the necessary changes in subscripts mentioned earlier. 3. Special Case of One Damper Winding on the Q-Axis The definitions have to be changed in this case. The absence of the g-winding can be expressed by assuming that: R = 0 (1-18) g and L = 0 0 (1-19) g Substitution of (1-18) and (1-19) into (1-11) yields the following result: T" = ^2. (1-20) qo RQ The substransient inductance L" then becomes: q ,M 2 W ^ v a"21) No transient quantities can be defined in this case. 103. APPENDIX 2 TRANSFORMATION OF THE EQUATIONS OF THE ELECTRIC PART The voltage equations of a synchronous generator in d,q,0-coordinates have the following, well known form: v = -[R]i - [L']i - [L ]±- i (2-1) -p — p -p p dt -p Application of the trapezoidal rule of integration yields: r | i (t) = 7 ^ i (t-At) - [L ] _ 1 v ( t ) - [ L n ] _ 1 [ R ] i ( t ) - [L ]_1[L«(t)]i (t) -A t —p A t —p p —p p —p p p —p ' [ L p ] ~ 1 v p ( t - A t ) - [ L p ] " 1 [ R ] i p ( t - A t ) - [ L p ] " 1 [ L p ( t - A t ) ] i p ( t - A t ) (2-2) Simple rearrangements allow rewriting of (2-2) into the following form: v (t) = [ R c o m p ] i (t) + hist (t-At) (2-3) -p -p p where the companion resistance matrix i s : [ Rcom P ] = _ f_2_• j + [ R ] + [ L ^ ( t ) ] ] ( 2 _ 4 ) and the past-history terms are: hist p(t-At) = [L ] - [R] - [L p(t-At)]]-i(t) - V p ( t - A t ) (2-5) Equation (2-3) has the form of a resistive companion model. This form i s preserved i f the equations are transformed from Park's coordinates to phase coordinates. The resulting equations can be written in the follow-ing way [2 ]: v(t) = [ R C O m p ]i(t) + hist(t-At) (2-6) — a,b,c where the subscript "a,b,c" denotes phase a,b,c-coordinates. The matrices [ R C O m p ] and [ R C O m p ] are not constant, and they have to be 3. f D y C recalculated at each time step. The amount of calculations i s , however, significantly smaller for the matrix [ R C O m p ] than for the matrix [R^° m P c], 104. i f the m a t r i x [L^] i s constant. This a f f e c t s the numerical e f f i c i e n c y of the s o l u t i o n . Because of t h i s , d,q,0-coordinates were chosen f o r the f i n a l a lgorithm. 105. APPENDIX 3 MULTIPHASE THEVENIN EQUIVALENT CIRCUIT OF A TRANSMISSION NETWORK I t i s assumed that a l l network parameters are l i n e a r . This i s the only assumption necessary to permit the c a l c u l a t i o n of the Thevenin equivalent c i r c u i t of a system shown s c h e m a t i c a l l y i n F i g . 3.1. N E T W O R K F i g . 3.1. Schematic r e p r e s e n t a t i o n of the network. The object i s to c a l c u l a t e the Thevenin e q u i v a l e n t c i r c u i t of the net-work seen from the terminals l a - 2 a , lb-2b, 1N-2N. The network i s described by the f o l l o w i n g nodal equation: [ G ] - v = i (3.1) where [G] i s a conductance m a t r i x created by a p p l i c a t i o n of the trap e -z o i d a l r u l e of i n t e g r a t i o n to d i f f e r e n t i a l equations of the network. The Thevenin equivalent c i r c u i t can be obtained as f o l l o w s : (1) S h o r t - c i r c u i t a l l voltage sources and cancel a l l current sources i n the network, (2) Connect a current source of +1.0 p.u. to t e r m i n a l l i and of -1.0 p.u. t o t e r m i n a l 2 i , and sol v e Eq. (3.1) f o r v. This 1a lb 2a IN 2b 2N 106. produces a column vector v which i s , in effect, the d i f f e r -ence of the l i - t h and 2i-th columns of [G] \ i.e. +1.0 in l i - t h component [G]v = £ except (3.2) -1.0 in 2-i-th component produces a vector: ^ T i = ^ l i - ^ 2 i ( 3 ' 3 ) The vector v m. i s the i-th column of the Thevenin resistive — T i • r_terminal, , _ matrix [R^ ] in Eq. (83). (3) Repeat step 2 for a l l other terminal pairs of interest, i = 2, ... N. This ends the calculation of the Thevenin re-sistive matrix [R^ e r n"" n a''"]. The open circuit voltages y^ 0(t) are calculated by the Transients Program as described in [45], 107. APPENDIX 4 RESISTIVE COMPANION MODEL The general property of i m p l i c i t i n t e g r a t i o n methods to create r e s i s t i v e companion models w i l l be demonstrated f o r the e l e c t r i c part of a synchronous generator. The procedure f o r other network components i s very s i m i l a r . The d i f f e r e n t i a l equations of a generator i n d,q,0-coordinates can be rewritten into the general form of: where: and [ C J = - [ L p ] 1 ( [ R p ] + [I/]) (4-2) [C ] = - [L p] 1 . (4-3) As i s w e l l known, the exact s o l u t i o n of (4-1) has the following general form [74]: i (t) = [ e [ C l ] A t ] - i (t-At) + f [ e [ C l ] ( t " T ) ] [ C ]v (T) dx (4-4) -P -P £ _ A t 2 -p Af t e r a p p l i c a t i o n of any i m p l i c i t i n t e g r a t i o n technique, (4-1) can be rewritten i n t o the following form: i (t) = [C l i (t-At) + [C ]v (t) + u (4-5) -P 3 ~~P h P — where _u represents the remaining part of the i n t e g r a l from (4-4) which contains only known "past h i s t o r y " values at t-At, t-2At, et c . The d e f i n i t i o n of the matrix [C 3] depends on the type of i m p l i c i t i n t e g r a -t i o n technique used i n the so l u t i o n of (4-1). For the trape z o i d a l r u l e of i n t e g r a t i o n , i t i s given as follows: [C 3] = ([I] - [C^r^ai] + ^| [ C J ) (4-6) 108. Simple rearrangements of (4-5) r e s u l t i n the f o l l o w i n g r e l a t i o n s h i p : v (t) = [C ] _ 1 i (t)' - [C r V ]i ( t - A t ) - [C ] _ 1 . u (4-7) —p it —p it 3 —p it — or i n a s h o r t e r form: V p ( t ) = [ C 5 ] i p ( t ) + h i s t ( t - A t ) (4-8) Equation (4-8) has the d e s i r e d form of a r e s i s t i v e companion model and i t has been obtained without s p e c i f y i n g the type of i m p l i c i t i n t e g r a t i o n method. Therefore, the r e s i s t i v e companion model can be c a l c u l a t e d f o r any system of equations of the form of (4-1) independent of the type of i m p l i c i t i n t e g r a t i o n technique used i n the s o l u t i o n of (4-1). The c a l c u l a t i o n of the m a t r i x [C 5] and the vec t o r h i s t ( t - A t ) i s very simple when the t r a p e z o i d a l r u l e of i n t e g r a t i o n i s a p p l i e d . APPENDIX 5 REDUCTION OF THE GENERATOR EQUATIONS ,red-| [ R r e d ] = [R ] - [R ][.R ] _ 1 [ R ] ss ss s r r r rs d The reduced r e s i s t i v e m a t r i x [R ] i s defined as f o l l o w s : ss (5-1) where-matrix [R ] has the f o l l o w i n g , w e l l known form ( i n d,q,0-s s coordinates) [ 8 ] , [ 9 ] : ss *33 a l l a 1 2 0 a 2 1 a 2 2 0 0 0 a 3 3 are fun c t i o n s of the (5-2) The elements a 1 ; L , a 2 2 , ar The elements a 1 2 and a 2 1 are f u n c t i o n s of both co and At, w i t h a l i n e a r dependence on co. Proof The m a t r i x [ R r r ] has the f o l l o w i n g form: [R ] = r r >1 I 0 I B, (5-3) where a l l the nonzero elements are funct i o n s of At only. The i n v e r s e m a t r i x [ R r r ] ^ i s simply [75] [R ] r r -1 -1 I B -1 (5-4) The matrices [R ] and [R ] are defined as f o l l o w s : rs s r [R 1 = r s 0 0 a 5 1 0 0 0 a 6 2 0 0 a 7 2 0 (5-5) with a l l the nonzero elements being functions of At only; [R ] = sr a ! 5 a i 6 a ! 7 a„, a„ a a lh 25 26 2 7 0 0 0 0 X5-6) with the elements a 1 1 +, a 1 5 , ,a 2 g, and a 2 7 being functions of At only and the elements a 1 6» a 1 7> a2h' a n d a 2 5 being functions of both co and At with linear dependence on co. Simple matrix multiplication yields: b. [R ] 1tR 1 = rr rs 11 12 32 hi (5-7) where a l l the nonzero elements are obviously functions of At only. Finally, the product [R ][R ] _ 1[R ] i s given as: SIT ITIT ITS [R ][R ] 1[R ] = sr rr rs > aik b l l + a ! 5 b 1 2 a 2 4 b l l + a 2 5 b 1 2 a i 6 b 3 2 + a ! 7 \ 2 a b + a b . 26 32 27 hi C l l c 1 2 0 — = C 2 1 C 2 2 0 (5-8) 0 0 o_ elements a i r , i 6 a i 7 ' a i h ' and a 2 5 depend (linearly) on co, the elements c 1 2 and c 2 1 in the resulting matrix w i l l have the same type of dependence. Substitution of (5-8) into (5-1) yields the desired result: ss which has the form of (81). l l l 21 l 1 2 22 33 (5-9) 111. ire d ' The m a t r i x [ R g s ] n a s t o be c a l c u l a t e d only once as the cons-tant of a s i m u l a t i o n run i f At does not change and i f s a t u r a t i o n i s not considered. The elements d^ 2 and d 2 1 depend l i n e a r l y on co , and t h e i r updating f o r changes i n co i s then obviously q u i t e simple. S i m i l a r procedures can be used f o r the four-phase e q u i v a l e n t c i r c u i t of the generator mentioned i n chapter 3.5. The r e s i s t i v e m a t r i x red [ R g g ] w i l l then be given as: [ R r e d ] -ss e l l 612 0 e21 e22 0 0 0 e33 e 0 0 •lk •21+ kk (5-10) where the elements e 1 1 9 e l l f, e 2 2 , e 3 3 , e ^ , and e ^ are f u n c t i o n s of At only. The elements e 1 2 , e 2 1 , and e 2 l t are f u n c t i o n s of both co and At w i t h l i n e a r dependence on co. I t can-also be shown that the f o l l o w i n g e n t r i e s of (5-10) are equal to the corresponding e n t r i e s of (5-9): (5-11) (5-12) (5-13) e i 2 = d12 e = d 22 22 e = d 33 33 112. APPENDIX 6 PRACTICAL CALCULATION OF THE MATRIX [R r e d] ss The Gauss-Jordan elimination process, which has been chosen, as the most efficient method, produces not only the reduced matrix [ R r e d ] , but also the distribution factor matrices [ D l = [R ]•[R 1 1 ss 1 2 sr rr and [D 1 = -[R ] - 1[R ]. The matrix [D; ] is needed for the calcula-2 1 rr rs 1 2 tion of the right-hand side of (80), and the matrix [I>21] i s necessary in the solution of the equations of the concealed terminals, once the retained variables have been found. A short description of the algorithm is followed by a flowchart of a computer program based on this algorithm. Consider a system of linear algebraic equations: [C]x = b (6-1) where: [C] = n x n matrix of coefficients, and x, b_ = vectors with n components. The objective is to arrive at a reduced system of equations for subset 1 in (6-2): (6-2) where [CJJ] and [ C 2 2 ] are matrices of order m x m and (n-m) x (n-m), respectively, and x^, b^ are vectors with m components andix^, b_2 vectors with (n-m) components. Elimination of variables x_2 from subset 1 results in the following system of equation: 113. red ' c c _ 1 | ^ 1 2 ^ 2 2 2 l 1*1 -C^C 2 2 2 1 b — 2 X - 2 — (6-3) Equation (6-3) has the d e s i r e d form a l l o w i n g the c a l c u l a t i o n of v a r i a b l e s ^ without' the need t o c a l c u l a t e the v a r i a b l e s x 2 , provided b 2 i s known. The transformation of the ma t r i x of (6-2) i n t o the matrix of (6-3) i s c a r r i e d out by exchanging - one at a time - the v a r i a b l e s x., b. f o r 3 3 m+1 <^ j j< n. This i s the Gauss-Jordan e l i m i n a t i o n process. For example, i f the v a r i a b l e s x., b. are to be exchanged ( v a r i a b l e s x.,,, b.,. , e t c . 3 3 J+l J+l have already been exchanged) and the c o e f f i c i e n t s from the l a s t exchange are C £ ^ ^ » then the j - t h row of (6-2) may be w r i t t e n down as f o l l o w s : C ( ° l d ) X + C ( ° l d ) X + 4- r ( ° l d K + r h + c • n x T c x • + .. . + c. . x. + c.,. . N b .,, + 31 1 J 2 2 J J j j ( j + D J+l Exchange of x^, b^ y i e l d s : c - J (old) (old) (old) c. 1 C. 1 - 1 C. 1 C. 1 + 1 3 J J j J (old) c. = c. . 3 33 where I f (6-5) i s r e w r i t t e n v i t h the c o e f f i c i e n t s ; + c. b = b . j n n j (6-4) ,(old) " . i n • c. 3 b - x. n 3 (6-5) (new) _ (new) , (new), , , (new), f r c. x + ... c... '.x. + c.. 'b. + ... + c. 'b = x. (6-6) J l 1 j ( j - l ) 3-1 33 3 n J then i t i s obvious that the f o l l o w i n g r e l a t i o n s h i p s h o l d f o r the e l i m i n -ated row: , „ (old) f o r j \ k (6-7a) (new) _ C j k "jk .(old) '33 and (new) ' j j ,(old) otherwise (6-7b) 114. Insertion of (6-6) into the remaining equations results in the following relationships: c(old) (new) = (old) _ (old) _jk • ^ _ g a ) lk ik i j (old) v J and c(old) (new) i j ^, o i _ \ c.. = — T ^ r r r otherwise (6-8b) 13 Aold) 33 A flow chart of a computer program executing the algorithm described above is presented on the next page. 115. START STOP A = 1.0/CU.J) K - 1 :i(K) = C(K,J) K = K + 1 = -C(J,K) x A DO 3 I = ], N 3 C(I,K) = C(l,K) + B x C1(I) DO 15 K 15 C(K,J) = = 1, N C1(K) x A C(J,J) J = .. = A -1 F i g . 6-1. Flow chart of computer program f o r matrix r e d u c t i o n . APPENDIX 7 CONVERSION OF THE OPEN CIRCUIT CHARACTERISTIC TO A FLUX-CURRENT CHARACTERISTIC The open-circuit characteristic V^j^g) = f ( i f ) 1 3 measured under balanced no-load conditions, i.e., i = i , = i = i _ = i _ = i = 0 (7-1) a b c D Q g The voltage equations (in phase coordinates) for these conditions are as follows: v = co Mj. sin (3 i , . (7-2) a r 1 r v b = co Mf sin 6 2 i f (7-3) v = lo M, sin B_ i £ (7-4) c f 3 r from which follows that sinusoidal changes i n flux are followed by sin-usoidal changes in voltages, since i f hysteresis and eddy current losses are neglected, as per assumption 5, chapter 5.4. Therefore, the RMS value of voltage difference between any two phases i s given as follows: co-M^-i. /3 V, , i_f_=«J± (7_6) (rms) 2^ /2 From (7-6) follows that the conversion of V^^g^ values to the instan-taneous flux values becomes a simple rescaling: = v ^ ! (7_7) CO The excitation current i ^ does not have to be converted. Consider the voltage equation for v , the only nonzero voltage under balanced no-load conditions: 117. V q = u;/2 M f # 1 f = * *d (7-8) Equation (7-8) can be visualized a system of coils placed in a rotating f i e l d create.d by a permanent magnet. This situation i s shown schemati-cally in Fig. 7-1. O Fig. 7-1. Schematic representation of an unloaded generator. The decrease in the value of (due to saturation) results in a de^ creased value of the flux linkage ty^, and therefore in a decreased vol-tage v . However, since no a.c.-components are present, there are no harmonics generated in the f i e l d distribution, i.e., the voltage v^ is s t i l l described by a linear equation of the form of (7-8). Therefore, the converted open circuit characteristicbhas the same form, as the original curve. APPENDIX 8 EXAMPLE FOR A MULTI-MASS SYSTEM The case presented here is a standard IEEE benchmark test case for the simulation of subsynchronous resonance phenomena [5]. The effects of a simultaneous three-phase short-circuit in a system as shown in Fig. 8.1 and 8.2 are simulated. GEN Rj =02 B R 0 - . 5 X, =i-14 Xj=j.5 Xc=:-j.371 INFINITE BUS X0=j-06 [Xf =j.04 Fig. 8.1. System diagram. Df RECTI'ON OF ROTATION Fig. 8.2. Model of the shaft system. The three-phase fault was applied at bus B at time t = 0 and then removed after time 0.075 sec, as soon as the current in each phase crossed zero. 119. The r e s u l t s presented here were obtained w i t h method I w i t h a time step of 100 usee. F i g . 8.3 shows the simulated electromagnetic torque of the generator, F i g . 8.4 shows the simulated mechanical speed of the generator r o t o r , and F i g . 8.5 shows the torque on the s h a f t between the generator r o t o r and the e x c i t e r . The i n c r e a s i n g o s c i l l a t i o n s of t h i s torque point out the reason f o r shaft damage which occurred i n a r e a l case from which the data of F i g . 8.1 and 8.2 were derived. - 4 . 0 F i g . 8.3. Simulated electromagnetic torque of the generator. 002 TIME (s) F i g . 8.4. Simulated mechanical speed of the generator r o t o r . F i g . 8.5. TIME(s) Simulated torque on the s h a f t between the generator r o t o r and the e x c i t e r . 120. The results presented here are practically indistinguishable from those in [5], which were obtained with a program developed in industry.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Synchronous generator models for the simulation of...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Synchronous generator models for the simulation of electromagnetic transients Brandwajn, Vladimir 1977
pdf
Page Metadata
Item Metadata
Title | Synchronous generator models for the simulation of electromagnetic transients |
Creator |
Brandwajn, Vladimir |
Date Issued | 1977 |
Description | Techniques for modelling of synchronous generators in the simulation of electromagnetic transients are described. First of all, an adequate mathematical model of the generator is established. It uses the conventional set of generator data only, which are readily available, but it is flexible enough to accommodate additional data, if and when such become available. The resulting differential equations of the generator are then transformed into linear algebraic equations, with a time varying coefficient matrix, by using the numerically stable trapezoidal rule of integration. These equations can be interfaced with the equations of an electromagnetic transients program in one of two ways: (a) Solve the equations of the generator simultaneously with the equations of a three-phase Thevenin equivalent circuit of the transmission network seen from the generator terminals. (b) Replace the generator model with a modified Thevenin equivalent circuit and solve the network equations with the generator treated as known voltage sources e[sup red][sub ph] (t-Δt) behind constant resistances [R [sup red][sub ph]]. After the network solution at each time step, the stator quantities are known and used to solve the equations for the rotor windings. These two methods cover, in principle, all possible interfacing techniques. They are not tied to the trapezoidal rule of integration, but can be used with any other implicit integration technique. The results obtained with these two techniques are practically identical. Interfacing by method (b), however, is more general since it does not require a Thevenin equivalent circuit of the network seen from the generator terminals. The numerical examples used in this thesis contain comparisons with field test results in order to verify the adequacy of the generator model as well as the correctness of the numerical procedures. A short discussion of nonlinear saturation effects is also presented. A method of including these effects into the model of the generator is then proposed. Typical applications of the developed numerical procedures include dynamic overvoltages, torsional vibrations of the turbine-generator shaft system, resynchronization of the generator after pole slipping and detailed assessment of generator damping terms in transient stability simulations. |
Subject |
Transients (Electricity) -- Mathematical models |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065485 |
URI | http://hdl.handle.net/2429/20635 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1977_A1 B73.pdf [ 5.59MB ]
- Metadata
- JSON: 831-1.0065485.json
- JSON-LD: 831-1.0065485-ld.json
- RDF/XML (Pretty): 831-1.0065485-rdf.xml
- RDF/JSON: 831-1.0065485-rdf.json
- Turtle: 831-1.0065485-turtle.txt
- N-Triples: 831-1.0065485-rdf-ntriples.txt
- Original Record: 831-1.0065485-source.json
- Full Text
- 831-1.0065485-fulltext.txt
- Citation
- 831-1.0065485.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065485/manifest