Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A synchronous generator model for the U.B.C. electromagnetic transients program Martinich, Terrence Gordon 1978

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1978_A7 M37.pdf [ 3.74MB ]
Metadata
JSON: 831-1.0065483.json
JSON-LD: 831-1.0065483-ld.json
RDF/XML (Pretty): 831-1.0065483-rdf.xml
RDF/JSON: 831-1.0065483-rdf.json
Turtle: 831-1.0065483-turtle.txt
N-Triples: 831-1.0065483-rdf-ntriples.txt
Original Record: 831-1.0065483-source.json
Full Text
831-1.0065483-fulltext.txt
Citation
831-1.0065483.ris

Full Text

A SYNCHRONOUS GENERATOR MODEL FOR THE U.B.C. ELECTROMAGNETIC TRANSIENTS PROGRAM by Terrence Gordon Martinich B.Sc. University of Victoria, 1972 B.A.Sc. University of British Columbia, 1976 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES in the Department of Electrical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA .May, 1978 c Terrence Gordon Martinich, 1978 In presenting th i s thes i s in pa r t i a l fu l f i lment of the requirements for an advanced degree at the Un ivers i ty of B r i t i s h Columbia, I agree that the L ibrary sha l l make it f ree ly ava i l ab le for reference and study. I fur ther agree that permission for extensive copying of th is thesis for scho lar ly purposes may be granted by the Head of my Department or by his representat ives. It is understood that copying or pub l i ca t ion of th is thes is f o r f i nanc ia l gain sha l l not be allowed without my writ ten permission. Department of p j f t r t i r i c a The Univers i ty of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T.1W5 Date May 8 , 13 78 i i ABSTRACT In r e c e n t y e a r s , prob lems have a r i s e n w i t h e l e c t r o - m e c h a n i c a l i n t e r a c t i o n in t u r b i n e - g e n e r a t o r u n i t s in power s y s tems . T h i s t o r s i o n a l i n t e r a c t i o n between the mechan i ca l s h a f t sy s tem and e l e c t r i c a l s y s tem o c c u r s in the subsynchronous resonance phenomenon, t h a t may a r i s e in s e r i e s compen-s a t e d t r a n s m i s s i o n l i n e s , and in f a u l t y s y n c h r o n i z a t i o n . S h a f t t o rques caused by the se d i s t u r b a n c e s may be s u f f i c i e n t l y s e v e r e as t o damage o r d e s t r o y the t u r b i n e - g e n e r a t o r s h a f t . T h e r e f o r e , to s tudy the se c a s e s , t h e r e was a need to r e p r e s e n t g e n e r a t o r s in d i g i t a l computer t r a n s i e n t s programs in more d e t a i l than as s i m p l e v o l t a g e s o u r c e s b e h i n d r e a c t a n c e s . A d d i t i o n of a synchronous g e n e r a t o r model to the U.B.C. T r a n s i e n t s Program i s d e s c r i b e d . The model o f the g e n e r a t o r r e p r e s e n t s both the e l e c -t r i c a l and s h a f t t o r s i o n a l dynamics and po s se s se s a l a r g e deg ree of g e n e r a l i t y . The model assumes c o n s t a n t e x c i t a t i o n and no magne t i c s a t u r a t i o n of the i r o n in the mach ine . I n t e r f a c e to the T r a n s i e n t s Program i s through a t h r e e -phase T h e v e n i n e q u i v a l e n t c i r c u i t o f the e x t e r n a l e l e c t r i c network as seen from the g e n e r a t o r t e r m i n a l s . T r a p e z o i d a l r u l e i m p l i c i t i n t e g r a t i o n i s ap -p l i e d to the machine e q u a t i o n s f o l l o w e d by use of Newton 's method to s o l v e the r e s u l t i n g n o n l i n e a r a l g e b r a i c e q u a t i o n s . T e s t cases a re p r e s e n t e d to v e r i f y the new s o l u t i o n t e c h n i q u e , f o l l o w e d by examples of p o s s i b l e a p p l i c a t i o n s of the synchronous machine m o d e l . In a s tudy of f a u l t y s y n c h r o n i z a t i o n s o f an un loaded g e n e r a t o r , c o n c l u s i o n s were in agreement w i t h the f i n d i n g s o f a m a n u f a c t u r e r . i i i TABLE OF CONTENTS Page ABSTRACT i i LIST OF ILLUSTRATIONS v ACKNOWLEDGEMENT vi 1 . INTRODUCTION. . 1 2 . SYNCHRONOUS MACHINE MODELLING . . k 2 . 1 Modell ing of E l e c t r i c a l E f fec t s k 2 . 2 Interface to the E l e c t r i c a l Network 11 2 . 3 Modell ing of Mechanical E f fec t s 12 3 . STEADY STATE - INITIALIZATION OF MACHINE EQUATIONS. . . . 16 3 . 1 E l e c t r i c a l Var iables 16 3 . 2 Mechanical System Var iables 18 4 . THE TIME STEP - NUMERICAL SOLUTION OF THE MACHINE EQUATIONS 20 4 . 1 Choice of Integration Method 20 4 . 2 So lut ion of the Machine Equations by Newton's Method. 20 4 . 3 Form of the Jacobian 23 5 . STEADY STATE AND TIME-STEP SOLUTION PROGRAMS . 3 1 5 . 1 Prel iminary Remarks About the Programs 31 5 . 2 Steady State - Subroutine "MSTEA" 33 5 . 3 Time Step - Subroutine , : IGEN" 3 ^ 5 . 4 Flow Diagram of So lut ion by Newton's Method 36 6 . TESTING OF THE SIMULATION METHOD 38 6 . 1 I n i t i a l Tests 38 6 . 2 S ingle Line-to-Ground Fault 39 6 . 3 IEEE Subsynchronous Resonance Benchmark Test Case . . 40 i v 7 . EXAMPLES OF PROGRAM APPLICATIONS J 43 7 . 1 Sustained Three-Phase Fault for an Unloaded Generator . . . . • 43 7 . 2 Faulty Synchronization of an Unloaded Generator . . . 43 8 . CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH. . . 48 REFERENCES 5 0 APPENDIX I 52 APPENDIX II 6 3 APPENDIX III 6 9 LIST OF ILLUSTRATIONS Fi gure Page 1 E l e c t r i c a l model of synchronous machine. 7 2 Terminal connection of the network to the machine. 12 3 Mechanical representat ion of a turbogenerator. 14 4 Model 1ing of mechanical system in the v i c i n i t y of the i - th rotat iona l mass. 14 5 Phasor diagram of the synchronous machine in steady s ta te . 17 6 Grouping of equations and independent va r i ab le s . 25 7 Nonzero s t ructure of the Jacobian matrix for an eight-mass seven-winding generator model (producing the largest matrix) 30 8 Flow diagram of so lu t ion by Newton's method. 37 9 Simple network conf i gura t ion . 38 .10 Network conf igurat ion fo r s ing le l ine- to-ground f a u l t . 39 11 E l e c t r i c network for the s imulat ion of subsynchronous resonance. 40 12 Mechanical model of the shaft system. 41 13 Subsynchronous resonance, phase c current (upper), E l e c t r o -magnetic torque on the rotor (center ) , and EXC-rotor shaft torque (lower). 42 14 System diagram for three-phase f a u l t . 44 15 Electromagnetic torque on rotor (upper), and LPB-rotor mechanical torque (lower). 44 16 System diagram for f au l ty synchronizat ion. 46 17 Torques due to f au l t y synchronizat ion (3= 120°) 47 18 Torques due to fau l ty synchronizat ion (3= 240°) 47 vi ACKNOWLEDGEMENT I would l i ke to extend my s incere thanks to my superv i sors , Dr. P. Fenwick and Dr. H. Dommel for the i r guidance, advice and unremitting encouragement and patience during my period of graduate study. I a l so wish to express my deep apprec iat ion to Dr. V. Brandwajn for his enthusiasm for this p ro jec t , his timely adv ice, and for supplying some of the test data used in this thes i s . The f i nanc i a l support of the National Research Council of Canada and the Un ivers i ty of B r i t i s h Columbia is g r a t e f u l l y acknowledged. 1. 1. INTRODUCTION This thes is serves three main purposes. The f i r s t purpose is to develop and then to add the electro-mechanical synchronous machine dynamics to the framework of the e x i s t i n g U.B.C. Electromagnetic Trans ients Program in the form of an accurate and reasonably e f f i c i e n t s imulat ion rout ine. The second is to show that resu l t s of the new program are in agreement with those of other methods. The f i n a l purpose of th i s thes is is to demonstrate that the program may be used to study ce r ta in problems of current importance. D i g i t a l computer programs s p e c i f i c a l l y designed for the s imulat ion of electromagnetic t rans ient phenomena in e l e c t r i c networks, have become a very powerful general-purpose tool in power system ana ly s i s . Some elements of the power system can be represented with a high degree of s oph i s t i c a t i on , Others are o v e r s i m p l i f i e d , e i t he r because more de ta i l ed representat ion is not warranted fo r the p a r t i c u l a r case study, o r , because ce r ta in nonl inear e f f ec t s are too complicated to model. Unt i l recent ly i t had been the prac t i ce to model a synchronous machine in electromagnetic t rans ients pro-grams as s inuso ida l voltage sources E | ^ a x s ' n (<JJt+<j>) behind impedances R + jXV , with R being the armature res i s tance and XV the d i r e c t - a x i s sub-s J d ' s • d t rans ient reactance. This approach neglects subtransient sa l iency and as-sumes constant rotor f luxes a f te r the t rans ient d isturbance. The simple machine model gives qu i te s a t i s f a c t o r y resu l t s during a few cyc les a f t e r i n -i t i a t i o n of a disturbance for ce r ta in studies invo lv ing fas t t rans ients [1]. Some types of t rans ients pe r s i s t for more than a few cyc le s . This is p a r t i c u l a r l y true in cases where, fo l lowing a disturbance such as a f a u l t , there is i n terac t ion between the mechanical turbine-generator system and the e l e c t r i c a l system. Examples include the torque amp l i f i c a t i on caused by sub-synchronous r e s o n a n c e , out o f phase s y n c h r o n i z a t i o n , and f a u l t y ^ s y n c h r o -n i z a t i o n . In the wake o f the two subsynchronous resonance, i n c i d e n t s o f 1970 and 1971, r e s u l t i n g i n s h a f t d e s t r u c t i o n a t the Mohave p l a n t i n s o u t h e r n Nevada [ 2 ] , u t i l i t i e s r e c o g n i z e d the need t o d e v e l o p more d e t a i l e d a n a l y t i c t o o l s t o s t u d y the t r a n s i e n t b e h a v i o r o f synchronous g e n e r a t o r s o v e r time spans l o n g e r than two o r t h r e e c y c l e s . In 1974, a j o i n t e f f o r t by S o u t h e r n C a l i f o r n i a E d i s o n and P a c i f i c Gas and E l e c t r i c i t y r e s u l t e d i n development o f a program c o n t a i n i n g b oth e l e c t r i c a l and s h a f t t o r s i o n a l dynamics o f synchronous machines. T h e i r machine model was implemented w i t h B o n n e v i l l e Power A d m i n i s t r a t i o n ' s e l e c t r o m a g n e t i c t r a n s i e n t s program (a much expanded v e r s i o n o f the U.B.C. program) but has r e c e n t l y been s u p e r s e d e d by a f a s t e r v e r s i o n [ 3 ] . T h i s machine model was n o t t h e o n l y one d e v e l o p e d d u r i n g t h i s p e r i o d [ 4 ] , The U.B.C. T r a n s i e n t s Program was p r e v i o u s l y m o d i f i e d [5] t o a c -commodate a synchronous g e n e r a t o r model, t r e a t e d as a t h r e e - p h a s e network branch w i t h n o n l i n e a r and t i m e - v a r y i n g parameters t h a t c o n n e c t s the g e n e r a t o r t e r m i n a l bus t o ground. The major r e s e a r c h e f f o r t was d i r e c t e d t o the i n -i t i a l i z a t i o n and p a r t i c u l a r l y the t i m e - s t e p s o l u t i o n , i m p l e m e n t a t i o n , and t e s t i n g o f the program. T h i s t h e s i s p roceeds w i t h development o f t h e synchronous machine e l e c t r i c a l and s h a f t t o r s i o n a l dynamics f o l l o w e d by a p p l i c a t i o n o f a new t e c h n i q u e f o r s o l v i n g the machine e q u a t i o n s [6] based on Newton's method. The a c t u a l F o r t r a n code i s shown f o r t h e s t e a d y s t a t e and t h e t i m e - s t e p s o l u t i o n s w h i c h was implemented under c o n t r o l o f t h e T r a n s i e n t s Program. T e s t cases a r e p r e s e n t e d t o e s t a b l i s h the r e l i a b i l i t y o f the new s o l u t i o n t e c h n i q u e and t o compare w i t h r e s u l t s o b t a i n e d by o t h e r methods. The t h e s i s then d e s c r i b e s examples where a p p l i c a t i o n o f the program i s used i n s t u d y i n g 3. problems such as the shaft torques caused by a three-phase f au l t and fau l ty synchronizat ion of an unloaded generator. F i n a l l y , conclusions and recom-mendations are given for fur ther research. 4. 2. SYNCHRONOUS MACHINE MODELLING 2.] Modell ing of E l e c t r i c a l E f fec t s A mathematical model of a physical device is often employed to obtain information about performance and c h a r a c t e r i s t i c s without having to test the actual device in the laboratory or in the f i e l d . The term "math-ematical model" implies an equation or system of equations ( l i near or non-l i nea r , d i f f e r e n t i a l or a lgebra ic) that descr ibes the bas ic physical behavior of the device. The usefulness and convenience of phys ica l device modell ing is apparent when one considers the d i f f i c u l t y and expense in removing from serv ice a large turbogenerator for test purposes. Any mathematical model of a physical device involves a ce r ta in degree of s i m p l i f i c a t i o n of r e a l i t y . The model s t ructure is derived from observations and physical understanding but there is always the p o s s i b i l i t y or impos s ib i l i t y of f i t t i n g parameters of the model to whatever data is a va i l ab l e . It is of no use to se lec t a highly deta i l ed and complicated model i f the parameters cannot be completely determined or estimated from test data. However, the degree of complexity of the model must be cons i s tent with i t s purpose. For example, one could not study the change in f i e l d cur-rent fo l lowing a terminal short c i r c u i t with a generator representat ion of s inusoida l voltage sources behind reactances. S i m i l a r l y , i t would be im-p rac t i c a l to employ a machine model representing such complicated t rans ient e f f ec t s as eddy currents in the rotor when the time period of study is so short that such phenomena play no ro le in the overa l l behavior of the physical system. To be of p r ac t i c a l a p p l i c a t i o n , a model of a synchronous generator should accommodate two kinds of rotor : 5. (1) Sa l ient pole rotor used in slow speed h y d r o - e l e c t r i c generating p lant s ; (2) s o l i d l y - f o r g e d steel round rotor in uniform a ir -gap turbogenerator un i t s . Sa l ient pole synchronous generators often have more than one pole pa i r ; the lower the head of water behind the dam, the larger the number of poles. A l so, these machines f requent ly have damper or "amort i sseur" c i r c u i t s . These damper c i r c u i t s , taking the form of brass or copper bars housed in the laminated pole faces, serve several useful purposes [7], inc luding the damp-ing of rotor o s c i l l a t i o n s fo l lowing disturbances such as short c i r c u i t s or switching. Uniform a i r -gap machines, on the other hand, are usual ly two pole and s lo t s are mi l led into the s o l i d round rotor body to support the f i e l d winding. With the above observat ions, the mathematical model of the syn-chronous machine w i l l be set fo r th with in the framework of the fo l lowing assumptions [5]: (1) A l l f l u x dens i t ie s are proport ional to the currents producing them so that there is no saturat ion or other nonl inear e f f e c t s . (2) The machine has a balanced, smooth s tator and the s ta tor windings are Y-connected with the neutral point grounded. (3) Eddy current paths in the rotor body of uniform a i r -gap machines or in the amortisseur c i r c u i t s of s a l i en t pole machines can be repre-sented by concentrated hypothet ica l windings, one in the d i r ec t axis (D) , the other in the quadrature axis (Q_) , and another hypo-the t i ca l winding (g) in the quadrature axis to represent deep f low-ing eddy currents in the case of round rotor machines. (4) Each s e l f and mutual inductance may be expressed as a sum of a 6. constant and a s inuso ida l funct ion of rotor angle 0(t) and/or 20(t). 0 is the angle between phase " a " and the d i r ec t axis of the f i e l d winding measured in the d i r e c t i o n of r o ta t i on . (5) Mutual inductances between armature and f i e l d and between armature and damper windings are equal. (6) During the time of the s imulat ion, f i e l d voltage is constant at i t s steady s tate value. (7) The d i r ec t i on of s ta tor current flow is out from the generator and into the external network. The cage - l i ke damper winding of s a l i en t pole machines cons i s t s of many c i r c u i t s carry ing d i f f e r e n t currents . Exact representat ion would re -quire a large number of concentrated hypothet ica l c o i l s but for most purposes i t is s u f f i c i e n t to represent the damper winding as one c o i l on the d i r ec t axis and another c o i l on the quadrature axis [8], For round rotor synchro-nous machines the eddy current paths in the iron can be considered as a cage winding having an i n f i n i t e number of bars. The model of the e l e c t r i c a l system cons i s t s of three i d e n t i c a l , symmetrical ly placed armature windings 1 a 1 , 1 b 1 , 1 c 1 and four unequal lumped rotor windings 1 f 1 , 1 g 1 , 1 D1 ,and 'Q_'. A l l rotor windings except the f i e l d 1 f are short c i r c u i t e d as they are not connected to voltage sources. A schematic representat ion of the synchronous machine is shown in F i g . 1. K i r c h ' h o f f ' s and F a r a d a y ' s laws may be app l i ed , y i e l d i n g one equation per c o i l in phase coordinates. 7= - [R ] T - ^ ? ( l) = -[R]T - {[L]T} (2) ijj is the vector of f l u x linkages [8]. The mathematical model o f the synchro-nous machine in phase coordinates is def ined by the seven f i r s t order d i f -7. F i g . 1. E l e c t r i c a l model of synchronous machine. 8. f e r e n t i a l equations represented by (l) or (2). This model has been tested b by others and resu l ts agree qui te well with f i e l d tests [k], [3] and should be adequate for most types of electromagnetic t rans ients studies [10]. frame of reference attached to the s tator and the rotor equations in a sep-arate frame moving with the rotor. The se l f - inductances of the s tator wind-ings together with the inter-winding mutual inductances, contain components that are functions of rotor angle. Owing to the rotor-ang le dependence of the inductance matrix, the d i f f e r e n t i a l equations given by (2) contain time-dependent parameters, making so lut ion rather awkward. A great s i m p l i f i c a t i o n resu l ts by using only one frame of reference and attaching this to the ro tor , thereby avoiding the time^varying inductance c o e f f i c i e n t s . The transformation of var iab les that accomplishes th is change is usual ly ca l l ed Park's t rans -formation [11]. It def ines a new set of s tator vol tages, cur rent s , and f luxes from a pro jec t ion of the actual winding var iab les on three axes; one along the d i r e c t axis of the f i e l d winding, c a l l ed the d i r e c t or " d " ax i s , a second along an axis perpendicular to the f i e l d winding, c a l l ed the quadra-ture or " q " ax i s , and the th i rd on a s tat ionary "0" ax i s . The transformation of s ta tor currents and voltages to the new coordinate system is defined by: v p = [T]v (3) In der iv ing (1) or (2), the s tator equations were expressed in a [ T ] i P where P 1 0 [T] (5) o; i COS0 cos (0-120°) cos(0+]20°) (6) /2 /2 /2 9 . which is s l i g h t l y d i f f e r e n t to Park's o r i g i na l transformation. Angle 6, as before, is the angle between phase " a " and the d i r e c t ax i s . The matrix [P] is orthogonal: [ P f 1 = [ P ] 1 and the transformation is power invar iant [12]. Rewriting (2) as v = - [R]T and applying the transformation gives [ 5 ] : v = - [L At - [R ]T - [L1 ] i P P J dt p p p p where •'! L c 0 0 fa fa o L c 0 0 0 ° ^ ° J " q 0 0 Lc 0 0 0 0 J dt ' = 10 L 0 q d 0 0 0 0 fa fa 0 0 0 0 42 q 42 q 0 0 0 0 L f M f 0 0 M f L D 0 0 0 0 M q 0 0 M q L 9 { [ P f L P ]> 0 0 V2 q \I2 q -fa fa 0 0 0 0 0 0 (7) (8) (9) (JO) (11) 10. or v d = -R i s d d . dt*d - u^q V q = -R * s q dt^q vo = -R s o " d . dt*0 v f = " R f f ' " d . dVf 0 = " R D D d , dt^D 0 = - R Q o_-~ d , dt^Q 0 = -R g g d tv g i nkages, i s def i ned by: = L d ! d + ^ ibt " f  +J L i q + q -.-\ *o - W 0 (12) (13) *f -/^f'd + Lf'f + Hf ! D =/fVd + Vf + L D ! D 0^.  =fh\  + VQ + Vq ij) =/|M i + M i + L i r g V 2 q q q Q g g Equations (12) together with (13) descr ibe the e l e c t r i c a l part of the synchro-nous machine model in Park's coordinates. In addi t ion to (12) and (13), i t is useful to develop the electromagnetic torque on the rotor . The power at the generator terminals is [13]: P = v , i , + v I d d q q Subst i tut ing from (12) for v . and v gives: 'd - ~ q d , . . d p e--"di * d + iqJt*q'+w(*di<,--+qid>-Rs<id + iq) .2 . .2, (1*0 (15) 2 Of the three groups of terms in (15), R s ( ' d + ' q ) represents the power loss in the s ta tor res i s tance; (i , + i is a measure of the time rate d dt d q at q of change of magnetic energy; and t*>(^ di q - i|» i j ) represents the power trans-11 3 ferred across the air gap from rotor to stator. The air gap torque, T r o t o r » can be defined as: Te = air-gap power  rotor mechanical rotor speed The mechanical angle, 8m, and rotor speed, as™, are related to the electric angle, 0, and frequency, w, of the network by: 0 = n 9 m (17-a) P co = n com (l7 -b) P where n^ is the number of pole pairs in the generator. Therefore, the elec-tromagnetic torque is: T 6 . = n (if/.i - ill i ,) (18) rotor p Td q - Tq d 2.2 Interface to the Electrical Network The seven first-order differential equations expressed by (12) are underdetermined in that they contain ten variables; seven currents and three voltages (v^ is a known constant). However, the stator components of v, i.e. V j , v^, and v^, can be eliminated. As the generator is usually con-nected to a larger electrical network, the network equations and synchronous machine equations are coupled and must be solved simultaneously. The inter-face to the network provides the constraint between stator voltages and cur-rents, allowing machine terminal voltages to be expressed in terms of stator currents. At the connection points of the machine to the external network, i.e. the three-phase bus, the currents injected into the network are the arm-ature currents i , i, , i while the bus voltages are the stator voltages v , a b c a Vb' V c ' a S s n o w n Fig. 2. These boundary conditions are a result of apply-ing Ki rchhoff's laws to the connection points. At any instant of time, there-fore, the effect of the network on the machine is identical to the effect of 12. a three-phase Thevenin equivalent c i r c u i t connected to the machine terminals, Or: (19) \M = C ^ N ^ + e N o ( t ) "N J N where the subscr ipt "N " denotes network quant i t ies and " o " denotes open ci rcui t quanti t i e s . j 1111 i 11111 GENERATOR 7 Thev "N 0 0 e c ELECTRIC NETWORK F ig . 2. Terminal connection of the network to the machine. A l l that remains is to transform (19) to d-q-0 coordinates, obtaining ex-pressions for v . , v , v„ in terms of i ,, i , i_. which can be subst i tuted into K d ' q ' 0 d ' q ' 0 (12). This procedure is shown irn Appendix 1. 2.3 Modell ing of Mechanical E f fec t s The synchronous machine mechanical system cons is ts of a ser ies of rotat ing components coupled together by shaf t s , as shown in F i g . 3- As with the e l e c t r i c a l part of the generator, the complexity of the mechanical model is con t ro l l ed by the amount of ava i l ab le data and the intended purpose of the model. Many of the damping constants, for example, are very d i f f i c u l t to obtain and are often not even'measured by manufacturers. It would be qui te 13. impract ical to develop a machine mechanical model w i th , for instance, so-ph i s t i ca ted nonviscous damping and nonl inear tors iona l spr ings. Even i f a l l the parameters could be measured, such add i t iona l complexity is unwarranted as i t does not contr ibute to the understanding of the bas ic system behavior. A l i near multimass-spring-damper system having up to e ight lumped masses was se lected as the most general mechanical representat ion of a syn-chronous machine. The actual number of lumped masses used in a given sim-u la t ion depends on the ind iv idua l case. This model accommodates a high pressure (HP) turb ine, an intermediate pressure (IP) turb ine, three low pressure stages (LPA, LPB, LPC) , the rotor , an e x c i t e r (EXC), and a feed-water pump (PUMP) or other piece of rotat ing machinery. However, app l i ca t i on of th i s model goes beyond the representat ion of a turbogenerator. It could equal ly well represent, for example, a motor d r i v ing one or more mechanical loads such as pumps, or a combination of a machine, turb ines , and mechanical loads on a shaf t . The fo l lowing assumptions are made in der iv ing the mathematical model of the mechanical system. (1) The system to be modelled contains up to e ight rotat iona l masses. (2) Mutual damping is considered unimportant and only speed self-damping is used. (3) The re l a t i ve angular de f lec t ions in the shaft sect ions are small enough so that the lumped masses may be connected by l inear tors iona l spr ings. (4) The period of s imulat ion is short compared to the time constant of This l im i t a t i on was imposed only because i t was f e l t that an eight-mass re-presentat ion is s u f f i c i e n t l y general to accommodate any actual mechanical system l i k e l y to be encountered. 14. HP IP LPA LRB LPC ROTOR . EXC PUMP F ig . 3- Mechanical representation of a turbogenerator. F i g . 4. Modell ing of mechanical system in the v i c i n i t y of the i - th rotat iona l mass. 15. the governor system so that input mechanical torques are constant. F i g . k shows the i - th rotat iona l mass and i t s v i c i n i t y . Newton's second law in rotat iona l form appl ied to each lumped mass y i e l d s a ser ies of second order ordinary d i f f e r e n t i a l equations. For the i - th mass: d 2 e. de. J . — — - + D . — - + K. . ,(e. - e. J + K. . . . (e . - e . . . * Tm T e i d 2 " D T M-1 i 1-1' M+1s i i+l) = T. - T. (20) . where J . = Moment of i n e r t i a of the i - t h rotat iona l mass. This could include the i n e r t i a of a sect ion of shaf t . Uni t s : kg-m 2 D. = C o e f f i c i e n t of speed s e l f (viscous) damping. Units : N-m-s/rad K.j= Spring constant of the shaft connecting masses i and j . Un i t s : N-m 6. = Angle of rotat ion of the i - th mass, measured in the d i r ec t i on of rotat ion in steady s tate from a f ixed reference. Uni t s : rad T T = The appl ied mechanical input torque on the i - th mass, tending to rotate i t in the +0 d i r e c t i o n . Units : N-m e . T. = The electromagnetic torque on the i - th mass, tending to rotate i t in the -0 di r ec t i on . Equation (20) may be expressed as two f i r s t order d i f f e r e n t i a l equations with G. and 6. chosen as s tate va r i ab le s . The f i n a l resu l t i s : i i de. i dt (21-a) - I - M T E d6. D. K.. 1 K T. - T dt i i i • 1 With " i " being summed over a l l the masses, equations (21-a,b) describe the mathematical model of the mechanical system. 16. 3 . STEADY STATE - INITIALIZATION OF MACHINE EQUATIONS •3.1 E l e c t r i c a l V a r i a b l e s The i n i t i a l i z a t i o n code used f o r the e l e c t r i c a l e q u a t i o n s was tha t used in imp lement ing Method I o f r e f . [ 5 ] , s i n c e both the new t i m e -s t e p s o l u t i o n t e c h n i q u e and t h a t o f Method I employ a T h e v e n i n e q u i v a l e n t o f the e x t e r n a l network and b a s i c a l l y s o l v e the same sys tem o f e q u a t i o n s . For the purposes o f i n i t i a l i z a t i o n , o n l y b a l a n c e d s teady s t a t e c o n d i t i o n s have been c o n s i d e r e d , where no n e g a t i v e and z e r o sequence components a re p r e s e n t . The a s sumpt ion i s c o r r e c t where a De l ta -Wye g e n e r a t o r t r a n s f o r m e r i s p r e s e n t ( z e r o sequence terms on the g e n e r a t o r s i d e o f the t r a n s f o r m e r a r e i d e n t i c a l l y ze ro ) and where the t r a n s m i s s i o n l i n e s a re c o n t i n u o u s l y t r a n s p o s e d [14] . With s t a t o r v o l t a g e s s p e c i f i e d by the u se r a t t ime t=0, the s t a t o r c u r r e n t s a r e found by the E l e c t r o m a g n e t i c T r a n s i e n t s Program fnom a pha so r s o l u t i o n o f the e n t i r e e l e c t r i c ne twork . These c u r r e n t s i , i , , and i cannot 3 D C be t r a n s f o r m e d to d - q - 0 components because the r o t o r a n g l e i s s t i l l unknown. When the machine i s o p e r a t i n g in s teady s t a t e , damper c u r r e n t s and z e r o sequence q u a n t i t i e s a re z e r o . In the s teady s t a t e : di , d i di , d X n d i _ di _ _ i = __s. = _ i = _P_ Q,= _g =0 (22) dt dt dt dt dt dt v ' From e q u a t i o n s (12) and (13) w i t h ( 2 2 ) : RD'D = 0 V Q = 0 Vg = 0 ( 2 3 ) o r i D = 0 i Q = 0 i g = 0 (2k) Vd " "Vd " "Vq ( 2 5 ) v = -R i + a i L . i . +0x^/1 f (26) q sq udd v 2 f r v f = - R f i f (27) 17. A phasor diagram may be constructed [15] from the known terminal condit ions of the machine (V , I , and the angle between them), as shown in F i g . 5. 3 3 The resul tant vector of V + R I + jX I l i e s on the quadrature axis as a s a q a demonstrated in ref . [5]. Angle 5, the pos i t ion of the quadrature ax i s , is found from: 6 = Tan lm(V + R I + jX I ) a s a g a Re(V + R I + jX I ) a s a J q a (28) If the reference axis corresponds to phase " a " of the i n f i n i t e system, then 6 is the load angle. Using the convention that the quadrature axis lags the DIRECT AXIS QUADRATURE AXIS REFERENCE AXIS F i g . 5- Phasor diagram of the synchronous machine in steady s ta te . 18. d i r e c t a x i s by 9 0 ° , the p o s i t i o n o f the r o t o r i s : e = -6 + \ (29) For b a l a n c e d s t e a d y - s t a t e o p e r a t i o n , d - q - 0 c o o r d i n a t e s a re r e l a t e d to phaso r q u a n t i t i e s by : i + j i , = /3Te" j 6 ( 3 0 - a ) q d v + j v , = / 3 \ /e~ J ' 6 ( 3 0 -b) q a where I and V a re RMS p o s i t i v e sequence phasors [5]. R e f e r r i n g t o the phasor d i a g ram: i d = / 3 | I | s i n (y -.6) (31 -a) i = /31 I I cos (y " <S) (31 -b) q v d = / J | v | s i n ( a - S) ( 3 2 -a ) v = /JI VI cos(a - 6) ( 3 2 -b ) q From e q u a t i o n ( 2 6 ) , the f i e l d c u r r e n t , i f , can be f o u n d . v + R i - X ,i , q s q d d , = -a /J (33) T h i s completes the i n i t i a l i z a t i o n o f a l l the independent e l e c t r i c a l v a r i a b l e s 3.2 Mechan i ca l System V a r i a b l e s In the s t e a d y s t a t e , a l l second d e r i v a t i v e terms o f 0 w i l l be z e r o and a l l the masses w i l l r o t a t e at the same a n g u l a r f r e q u e n c y , u i m , g i ven by (I7"b). E q u a t i o n (20) may be r e w r i t t e n in the fo rm: K.. . ( 0 . - 0. J + K, . . , ( 9 . - 6 . x 1 ) = TT - T? - D . o j m , i = l , . . . , n / , M i i -1 i i+1 i i+l i i i { Q (34) nt»c> Upon add ing a l g e b r a i c a l l y the n e q u a t i o n s e x p r e s s e d by (34), (where n i s the number o f lumped masses on the s h a f t ) , a l l terms in K. d i s a p p e a r , l e a v i n g : 6 = E ( T m - T e - D.u)m) (35) 1=1 1 1 1 R e a r r a n g i n g g i v e s : I T m = S ( T e + D.u>m) (36) i=1 ' i=1 ' 1 T h i s g i v e s the t o t a l i npu t mechan i ca l t o r q u e d e v e l o p e d by the t u r b i n e s to b a l a n c e the o p p o s i n g e l e c t r o m a g n e t i c and damping to rques in the s teady s t a t e . e e The e l e c t r i c a l t o r q u e , T . , i s nonzero o n l y f o r the r o t o r , where T i s i ' r o t o r de te rm ined by e q u a t i o n s (13) and ( 1 8 ) . The mechan i ca l t o r q u e s , T m , a re nonzero o n l y f o r the t u r b i n e s . The t o t a l mechan i ca l i npu t t o rque i s c a l c u l a t e d from ( 3 6 ) . Knowing the p e r -centage o f the t o t a l mechan i ca l t o rque d e v e l o p e d by each t u r b i n e , the i n d i -v i d u a l t u r b i n e t o rques can be a s s i g n e d . Hav ing a l r e a d y d e t e r m i n e d the e l e c t r i c a n g l e o f the r o t o r , the p h y s i c a l p o s i t i o n o f t h i s ang le may be found from e q u a t i o n ( I 7 ~ a ) . W r i t i n g e q u a t i o n (3*0 f o r a l l the masses e x c e p t the r o t o r , the r e s u l t i n g s e t o f (n-1) l i n e a r a l g e b r a i c e q u a t i o n s i s e a s i l y s o l v e d f o r the rema in ing 6. . T h i s com-p l e t e s the i n i t i a l i z a t i o n p r o c e d u r e . 20. 4. THE TIME STEP - NUMERICAL SOLUTION OF THE MACHINE EQUATIONS 4.1 Choice of Integration Method The coupled f i r s t - o r d e r d i f f e r e n t i a l equations (12) and ( 2 1 ) , re -presenting the complete synchronous machine model, are solved simultaneously with the rest of the external e l e c t r i c network. A large number of numerical integrat ion schemes have been devised for so lv ing s tate equations but not a l l of them are su i t ab le for so lv ing the machine equations. T rapezo ida l -rule i m p l i c i t integrat ion was adopted for the fo l lowing reasons: (1) Trapezoidal integrat ion is numerical ly s tab le for any s teps i ze , At , for a set of s tab le d i f f e r e n t i a l equations [ 1 6 ] . (2) This second order method is s e l f - s t a r t i n g , requir ing no past h i s tory terms except from the immediately preceeding time step. (3) The choice of At is not a f fec ted by the smal lest time constant in the network; i t is r e s t r i c t e d only by des ired accuracy of the s o l u -t i on . (4) There seems to be no d e f i n i t e evidence that higher order i m p l i c i t integrat ion methods lead to numerical super i o r i t y [ 6 ] . (5) Use of the trapezoidal rule of integrat ion s i m p l i f i e s i n te r fac ing with the Electromagnetic Transients Program, which uses the same s teps ize and the same integrat ion technique. (6) This method is simple to program and is the s implest of a l l i m p l i c i t integrat ion methods [ 1 8 ] . 4.2 Solut ion of the Machine Equations by Newton's Method For the system of s tate var iab le equations: x = f (x , t ) (37) 21. a p p l i c a t i o n o f the t r a p e z o i d a l r u l e r e s u l t s in the a p p r o x i m a t i o n : x ( t ) = x ( t - A t ) + £1 { f ( x ( t ) , t ) + f ( x ( t - A t ) , t - A t ) } } (38) U s ing t h i s a p p r o x i m a t i o n , the d i f f e r e n t i a l e q u a t i o n s f o r the e l e c t r i c a l sy s tem become the a l g e b r a i c e q u a t i o n s : ^d = * d + {'Vd 'VV 6 0 " *d " h^d^c? ~ ("V^ q* \1> = $ + ~ ' {-v ( i ,, i , i - , 0 ) - v - R ( i +i ) + (wip +u$ ) r q y q 2 q v d ' q ' 0 ' q s q q r q T q ~ At ~ ^0 = ^0 + ~T { " V O ( i d ' ! q ' i O ' 0 ) " V 0 " R s ( i 0 + i 0 ) } *f = *f + { " ( y f + V " R i T f + ^ f > > (39) *D = + ^ { - R D { 1 D + r D ) } * = If + —Y (-R ( i + i ) } 9 9 2 g v g g ' where 0 i s the e l e c t r i c ang l e and co the a n g u l a r f r e q u e n c y o f the network , r e l a t e d to the r o t o r mechan i ca l a n g l e and speed by ( I 7 _ a , b ) . The m e c h a n i c a l e q u a t i o n s f o r the i - t h mass become: At_ 2 At . + =j ( e . + e.) (^ o) i " 9 i + I J - . { - D i 6 i + K i i - I e i - ! " ( Kii"l + Kii + 1 ) 6 i + K i i + 1 6 i - 1 i + T? - "1 - V i + K i i - i ^ - i - ( K i i - i + " W V Kii +iVi + T m - T 6 } (41) The s u b s c r i p t , " i V , runs f rom 1 to n, where n i s the number o f lumped masses on the t u r b i n e - g e n e r a t o r s h a f t . For c o n v e n i e n c e , a l l v a r i a b l e s at the o l d t ime t - A t , have been w r i t t e n w i t h " ~ " o v e r them, those a t the new t ime t have been w r i t t e n w i t h no argument " ( t ) " . E q u a t i o n s (39),(40) , and (41) a r e n o n l i n e a r a l g e b r a i c e q u a t i o n s 22. and an i t e r a t i v e so lu t ion is required. Very often Newton's method is employed in such cases because of (1) ve ry good e x p e r i e n c e w i t h i t in network power f l o w s t u d i e s [19] » and, (2) i t s convergence is quadrat ic and therefore fas t [20]. To i11ustrate the use of Newton's method, assume a so lut ion x ^ , . . . , x is required for the set of "m" nonl inear equations: g 1 ( x l ' - - - x m ) = 0 g 2 ( x T , . . . x m ) - 0 Newton's method proceeds as fol lows (i) (42) (1) Set i t e r a t i on count, i , to 0 and make an i n i t i a l guess, x , to the so lut ion vector, (2) Using x ^ , evaluate g { x ^ } and the Jacobian matrix [ J ^ ] , defined as: [J] = 3x, 3x, 3x m 9g 3g 3g, 3x, 3x, 3x x = x ^ \43) (44) 1 2 m (3) Solve the set of l i nea r equations for the vector Ax: [ J ( i ) ] A^= - g ( x ( i ) ) (4) A better approximation to the so lut ion is then found from: x ( i + 1 ) = x ( i ) + Ax (45) (5) If any one of |Ax^|, . . . , | A | are larger than a predetermined s, then increment the i t e r a t i o n count by 1 and return to step ( 2 ) ; otherwise the method has converged. 2 3 -Newton's method is appl ied to the machine equations to solve for the currents and mechanical angles at each new time step. Depending on the s t a r t i n g po int , the procedure may converge on a number of d i f f e r e n t so lu t i ons , only one being the des ired so lu t i on . In t rans ients s imulat ions, the s tep-s i z e , At , is normally less than one mi l l i second and under these circumstances extraneous so lut ions w i l l not be generated by Newton's method. 4.3 Form of the Jacobian Matrix Solut ion by Newton's method requires the synchronous machine equations to be wr i t ten in the constra int form of equation (42). Let g™ and g 6 be defined as the vectors of mechanical and e l e c t r i c a l const ra int equations, respect ively. Taken together, they form the r ight hand s ide of equation (44). Beginning with the mechanical system, equation (40) g ives: Subst i tut ing (46) into (4l) and rearranging terms gives one a lgebra ic constra int equation for each mass. For the i - t h mass: 9? = (-TT + ~~T—) ® • " o~T~~ ^ K. . ,6. . - (K. . ] +K.. J . 1 )6. + K .^ ,0 .^ , i At J j i 2J; - - i i - l i - l u -1 i i + l i i i + l i + l The e l e c t r i c a l equations contain the network frequency, u>. This frequency can be expressed initerms of the rotor angular frequency from (40) with (17-b) as: u = n 8 = n tvr- 9 - ( § + 0 )'] (48) p r p L At r v r At r' v ' where the subscr ipt ' V refers to the rotor and n is the number of pole 2k. pairs in the synchronous machine. Using equations (12) and ( 4 8 ) , the e l e c t r i c a l cons t ra in t equations are def ined by: . A - * d ( i d . V D ) - ^ + - T - t v d ( , d ' , q ' , o ' e r ) + V R s ( i d + y } *Vq ( , q V g ^ 9 " - " ^ - " ^ q V V V ' V = ° 4 - W W " *q + ^ f " { v q ( i d ' i q ' i o ' 6 r ) + V R s ( i q + i q ) } - n p * d ( , d ' , P i D ) , ( e r " 5 r ) + - T - " [ * d ( 1 d ' 1 f ' D>-*d ] = 0 93 " • o ( , o ) - * o + { v o ( i d » i q . i 0 ) + v 0 ^ s ( i 0 ^ 0 ) > - 0 g^ = * f ( i d , i f , i D ) " * f + -f-<VVRf(if+'f)} = 0 95 - W ' P ' D ^ D + " T : { R P < I D + ' D ) } = ° 96 - V v W _ i G + ^ V W } = ° y 7 V q 0., g g 2 . g g g Before applying Newton's method to the a lgebra ized machine equations, i t is advantageous to arrange them in the fo l lowing groups (F ig . 6 ) . Except for the rotor , Group 1 of the vector " g " contains a l l the mechanical equa-tions beginning with the H.P. turbine and ending with the e x c i t e r and the feed water pump. Group 2 contains a l l the e l e c t r i c a l equations s t a r t i n g with the d-q-<j> voltage equations followed by the voltage equations for the rotor c o i l s and, f i n a l l y , the mechanical equation for the rotor . The elements of the vector of independent va r i ab l e s , " x " , are arranged as shown. (49 - a ) (49 - b ) (49 - c ) (49 -d ) (49 - e ) ( 4 9 - f ) (49 - g ) 25. CONSTRAINT VECTOR INDEPENDENT VARIABLES 9 = m 9l m r-1 m r+1 m 'r+2 e 9l e GROUP 1 GROUP 2 '1 \ 3 r-1 3r+l 3r+2 i F i g . 6. Grouping of equations and independent va r i ab le s , " r " ind icates rotor quan t i t i e s . The Jacobian elements corresponding to Group 1 are: — = - — • K (50) 30. . 2J. Ii—l ^ ' i - l I 99? . , D, i - w r ~ ar + i. + I T . - ^ M - I ^ I W ! ( 5 , ) I l l 3q m (i + l i Pa r t i a l der iva t i ves with respect to a l l other independent var iab les are zero. 2 6 . Taking p a r t i a l d e r i v a t i v e s , the Jacobian elements for Group 2 are as fo l lows: I R S + ] (53) l p e r -( 5 ^ t + n p e r ) ] . L q (5*) - L + A t 9 i d L d ' 2 3g* At 9 v d ai q 2 3i q 3v At 3i o 2 3i 1 e .39], e 9g 1 3 i f 9 i D 3g| ' ag* 9 V 17 ai g (55) \ M f (56) 9 gt 9 gt At 9 v d l e 7 = n P 1 F - = n P [ — 1 e — + V ( 5 8 ) 8 9 2 " A t 3 " V q " " [n 6_ -( & + n i ) ] . L , (59) 8 i d 2 3 i d p r v w 2 p r' 39o . 3 v — ^ - L + - ^ [ R + ^ T 1 ] (60) ai q 2 s ai q q 9 g 2 _ At % ai 2 ai o o e e 3 i f 3 i n w 2 .V 2 " f 992 992 3i_ ai v 2 q Q g (61) " 2 - 9 g 2 _ ^^[J_, m - y i M n (63) _ 1 = r At _ q _ _ j ^ ( 6 i | ) 3G p L 2 3© V d J r • 9 93 At 3 V0 2 393 At 9 v o 3i q 2 3i q >«;' 1; + ^ 3 i o (65) (66) 9 v n + J T ] ( 6 7 ) (The pa r t i a l der iva t i ves of v d > v , v Q with respect to i d , i q , i Q are c a l -culated in Appendix I). U + *lR' (68) 3 i f " 2 " f 3 9 c 5 - U + ^ R „ (69) 3 i D D 2 D 3 i d 3 i d VZ f 5 = Jl M , (70) 994 3 g | 3 i D 3 i f f dH . . At = M r (71) 3 9 ' - L • ^ R (73) 3i g 2 g 9 39^ 39y 3i 3 i v 2 q q q 9 9 6 9 9 6 a i 3 i n — = M (75) Ror the rotor : 2 8 . K - V l j V V V - K r r - l ( 6 r - l + 9 r - l ) + ( K r r - 1 + K r r + 1 >' < e r + 6 r> + K At I r rm , ;ra rr+1 ( e r + 1 + e r + 1 ) - ( T ; + T ; ) + ( T ^ T p } =o (76) Electromagnetic torque, , is known from (18) S p " § r + U { D r ( e r + § r > " K r r - 1 < e r - l ^ r - l > + < K r r - 1 + K r r + 1 > * W r + K r r + 1 ( e r + 1 + e r + 1 ) - (T^ + ^ ) + . . n p [ * d ( i d , i 1 ! , i D ) - i q + * d - i q ] - W ' q ' W ' d + V ' d ] } = ° (77) Subst i tut ing f o r 6 r from (41) g ives: ^ - + 17/ - K r r - l 6 r - l - + ( K r r - 1 + K r r + l > ' 9 r " " W r + l " ^ 2 D n -+ W d ^ f ' ^ ' q " * q ( , q » , Q ' , g ) ' , d ] } " [ A t + T . 1 ^ " " r + f j { ^ r r - l V l + ( K r r - l + K r r + l > ' K " K r r + l i r + l " ^ + n - i|i • ! . ] } = 0 (78) p d q q d Rrom ( 7 8 ) , the las t row of the Jacobian matrix i s : 41 K , (79) m 3 g r 9 6 r _ 1 2 J r ' V r - l m = - ^ K (80) ae r + 1 2 J r K r r + 1 a 9 ; [ L , _ ^ ] (81) 3 i d 2 J r j p L d q vq m ^r_= [ L j ] (82) 91 2J p yd q d q r r m ^ m 9 i f 3 i D 2 J r p LV2 f q 29. 3 9 r 3 g r A t -3 i4 " 9 i g " 2 J r " P r r r A l l other p a r t i a l s are zero. The nonzero s t ructure of the Jacobian, [ J ] , is shown in F i g . 7-The reason for arranging the equations as in F i g . 6 is now more ev ident, the Jacobian becomes pa r t i t i oned into two d i s t i n c t parts . The f i r s t is the time-independent par t , derived from the mechanical equations without the rotor. This part comprises 12% of the matrix in the case of an eight-mass seven-winding machine model, that i s , when the matrix is as large as poss ib le . The time-dependent block on the lower r ight comprises only 28% of the matrix for th is case. As the number of lumped masses in the model decreases, so does the s i ze of the time-independent part of [J] u n t i l , for the l im i t i n g case of a one- ( i n f i n i t e ) mass system, there is no time-independent part at al 1 . For s imulat ion purposes, some computer time is saved by exp lo i t i n g the s t ructure of the Jacobian. The time-independent part , obv ious ly , need not be updated during success ive i te ra t ions or from one time step to the next. A l so , by wr i t i n g the equation for the rotor , l a s t , the number of e l iminat ions is reduced during LU t r i a n g u l a r i z a t i o n of the matrix p r i o r to the back sub-s t i t u t i o n process for so lv ing the AX.. CD: < < T tn (O < <_ t— c _ c _ <_ . c _ O O v n U J —* < < < < CD : CD CD : CD : CD CD : a>' -ti O Xl Q_ oo U l - f U J ro —* II II II II ll II II II II ll ll 3' .7+ O O 0--:-s (D — QJ (D 3 - i ^ O T J CD -I — O l O D_ zy ~\ C r t C n i o — 3 r t 3 QJ C to in - i in (D r t 3" in O (D (T> - t < — ' (D r t Q) 3 3 -I I (D O " J N in in 3 oj r t a n — o 3 3 cr QJ I Q — r t QJ -I I Q 3 — CD X 3 3 CO QJ • -1 r t QJ -1 r t — O X r - X X X X X m. o P X X X o m -D ^> X X X m z o m uP X X X z H ~P X X X coCD xl X X X X X X X X X X X X X H rn c — • O X X X o m X X X X X ~o rn z o X X X X X rn z H X X X X X L O X X X X X X X X X X X 1 • TIME DEPENDENT U J o TIME INDEPENDENT 3 1 . 5. STEADY STATE AND TIME-STEP SOLUTION PROGRAMS 5.1 Prel iminary Remarks About the Programs The fo l lowing information must be suppl ied by the user i f a synchro-nous generator model is to be used in the Transients Program. (1) the e l e c t r i c a l frequency. (2) the peak l i n e - t o - n e u t r a l voltage at the generator terminals and the angles of the three phasor voltages with respect to the reference. 7*3 f3 TjM r, L,. Lr,, L , , L . Ln> and L r t in h. 2 f f D' q'v2q g Q 0 (4) the f i v e res istances R , R £, R^, R_, and R in ft. s f D Q g (5) the number of poles in the generator, the number of lumped masses on the sha f t , and the mass number corresponding to the rotor. (6) i ne r t i a s of the lumped masses, in kg-m 2. (7) spring constants for the shaft sec t ions , in N-m. (8) damping constants, in N-m-s/rad. (9) the percentage of the tota l mechanical torque contr ibuted by each turbine stage. In addi t ion to the above, small zero sequence and pos i t i ve sequence res istances must be s p e c i f i e d so that, in the steady state c a l c u l a t i o n s , the Transients Program can convert the generator to a Norton equiva lent c i r c u i t . A small but f i n i t e res is tance (»0 .001 f t ) is required across each current source. This program converts the two resistances to s e l f and mutual quant i t ie s which are subsequently included with the s ta to r res istances of the machine during steady state i n i t i a l i z a t i o n . Even the smal lest coupling res is tance permitted may be large compared to the actual armature res i s tance. To overcome th i s problem, the zero and pos i t i ve sequence res istances should be set equal to each other, for instance, 0.001ft. By spec i fy ing an armature res is tance equal 32. to the true value minus 0.001ft, the f i c t i t i o u s res is tance disappears when absorbed into the s t a to r . The program assumes that, in l a b e l l i n g the masses on the shaf t , the high pressure turbine appears f i r s t , fol lowed by the intermediate pressure stage and so on down the shaf t . It may be that the user does not wish to model the mechanical system but requires the e l e c t r i c a l desc r ip t ion of the synchronous machine. By spec i fy ing a two-pole, one-mass system with the rotor corresponding to mass number 1, the timesstep program solves for a machine with an a r b i t r a r i l y large rotor mass having constant rotat iona l frequency. Quite o f ten , the manufacturer's data of the e l e c t r i c a l parameters of the machine inc lude: X,, X ' , X'!, T' T 1' , X , X1 , X", T' , T" , and R c. d d' d dO dO' q ' q q qO qO f The funct iona l re l a t ionsh ip between these constants and the desired set of res istances and s e l f and mutual inductances is par t ly nonl inear. The param-eters required by the program may be extracted from th i s data by applying Newton's method, described in ref. [ 5 ] . Sometimes T' T'J, T ' , and T" are ' 1 1 dl d ' q ' q given instead of the open - c i r cu i t values. The desired time constants for the d i rec t axis can be obtained with [21]: X X T io " 7 1 T i <86> d K 'dO Y 1 1 'd (87) X d Extension to the quadrature axis is obvious. If L is set equal to zero, there is no " g " winding on the quad-rature axis and the program solves for a s ix-winding e l e c t r i c a l machine mode 1. 33. 5.2 Steady State - Subroutine "MSTEA" The fo l lowing steps ind icate the bas ic flow of events during the steady state i n i t i a l i z a t i o n of machine currents and mechanical angles by subroutine "MSTEA11. I I n i t i a l i z e E l e c t r i c a l Equations (1) Ca lcu late real and imaginary parts of p o s i t i v e , negat ive, and zero sequence currents from real and imaginary parts of the terminal currents given by the Transients Program. (2) The magnitude of negative and zero sequence currents indicates the degree of unbalance in the system. Both should be very small com-pared to the po s i t i ve sequence current . (3) Ca lcu late the real and imaginary parts of V + R I + jX I s a q a (k) Find the angle of the quadrature ax i s , <5: & = Arctan lm(V + R I + jX I ) 1 a s a J q a Re(V + R I + jX I ) a s a q a (5) Rotor angle £ = 6 j (6) Compute the magnitude of V and i t s angle, a. V is the RMS l i n e - t o -3 3 ground voltage for phase " a " . (7) Solve for voltages: v d = ^ l V a l 5 i n ' a " 6) v = /3|V | cos (a - 6) C] 3 (8) Compute magnitude of I and i t s angle, y. 3 (9) Solve for currents : i d = / J | 1 a | sIn(y " <5) i = fi\ l a | c o s ( Y - <S) (10).Calculate the f i e l d current , i ^ . 34, v + R i - ( i iL.i , q s q d d A l l other currents are zero. (11) Form Park's inductance matrix, f -p l * (12) Transform voltages and currents to phase coordinates fo r p r in tou t . (13) I n i t i a l i z e the vector of f lux l inkages, ip = [L^ ] i . II I n i t i a l i z e Mechanical Equations (1) Ca lcu late electromagnetic torque on the rotor plus the sum of a l l speed self-damping torques. (2) A l l oca te steady state input torque to each turb ine. The sum of the turbine torques equals the torque ca lcu la ted in Il-(1). (3) Form the c o e f f i c i e n t matrix, [K] , of the mechanical equations, omitt ing the rotor equation. (4) Form the r ight hand s ide vector g. = T1? - T? - D.cjm, again omitt ing the rotor. (5) Solve the set of l i near equations [K]6 = g for the 6.. These are the i n i t i a l mechanical angles at time t = 0. (6) Ca lcu late the mechanical h i s tory vectors : HISTl( i ) = 2 u>m + (^- + ji-)8. i HIST2 ( i ) = | T t -K . . i e - 1 + ( K - - 1 + K . . . J e . - •<.._,, 9. - T T + T?} 2J. i-1 M-1 1 i+ l 1 11 + I i + 1 1 1 1 (7) Return to the Transients Program. 5.3 Time Step - Subroutine "GEN" The t ime-step so lut ion of the synchronous machine equations cons ists of the fo l lowing steps: (1) P red ic t new mechanical angles 8. = B. + 9.At and new currents i = i . 35. ( " „ " means from the previous time s tep ) . In steady state operat ion, these pred ict ions are co r rec t . (2) Modify s ta tor res i s tance vector to absorb the small external coup l -ing res i s tances . (3) Compute the time-independent constants appearing in the Jacobian matrix, "JACOB", and the vector " G " . Terms in quotes appear in the actual Fortran code.in Appendix III. (k) Compute the time-independent elements (mechanical terms) of the Jacobian. (5) Ca lcu late the rotor e l e c t r i c a l angle and i t s angular frequency from past h i s tory and the predicted new mechanical angle. (6) Compute the time-dependent factors appearing in the Jacobian and the vector " G " . (7) Ca lcu late terminal voltage and f lux l inkage vectors corresponding to the present estimate of currents and rotor angle. (8) F i l l in the elements of the time-dependent part of the Jacobian. (9) Ca lcu late the vectors of the mechanical part that w i l l become h i s tory terms in the next time step. (10) Ca lcu late the r ight hand s ide vector "G " of mechanical and e l e c t r i c a l equations. (11) Solve for the correct ion terms, AX, in Newton's method from the set of l i near equations: [j] .AX = - i (12) Update the currents and mechanical angles: i j k L . f - 0 + A i . eK + A0J. j j j 36. ( superscr ipt "k" indicates i t e r a t i on number; k=1 ,2 , . . . 10 ) (13) If more than ten i te ra t ions have been performed, e x i t from the pro-gram with a warning message. The procedure is not converging. (14) If any of |AX.| are larger than predefined tolerance va r i ab l e , eps i lon (=10 9 ) , go to (5) and i te ra te once more. (15) Since convergence is successful at this po int , update terminal voltages, f l ux l inkages, and mechanical h i s tory vectors . (16) Store a l l quant i t ies needed for the next time step so lu t i on . (17) Convert d-q-0 currents to the phase domain to be in jected into the external network. (18) Return to the Transients Program. 5.4 Flow Diagram of So lut ion by Newton's Method Under the control of the electromagnetic Trans ients Program, the time step program solves the synchronous machine equations during each time step of the s imulat ion. F i g . 8 shows the flow diagram of the so lu t ion algorithm by Newton's method, with the in te r fac ing between the trans ients and machine programs in phase coordinates. START t = 0 t = t + At TRANSIENTS PROGRAM SOLUTION WITHOUT GENERATORS PREDICT CURRENTS AND MECHANICAL ANGLES USING THEVENIN EQUIVALENT OF NETWORK PROVIDED BY TRANSIENTS PROGRAM, CALCULATE CORRECTIONS BY NEWTON'S METHOD F ig . 8. Flow diagram of so lu t ion by Newton's method. 38. 6. TESTING OF THE SIMULATION METHOD The main purpose of running test cases was not to ve r i f y the syn-chronous machine model, s ince this has already been done by others [5], but rather, to confirm the p a r t i c u l a r s imulat ion method chosen. The tes t ing was done in various stages. 6.1 I n i t i a l Tests the Transients Program could be replaced by a Thevenin equivalent c i r c u i t . The network, shown in F i g . 9, cons isted of phase res istances of 10.fi con-nected to s inuso ida l voltage sources with amplitudes of 5000. vo l t s . Shaft tors iona l e f f ec t s were neglected by modelling the mechanical system as one a r b i t r a r i l y large rotat ing mass. The e l e c t r i c a l parameters of the machine were those used in Example k of ref . [5] and peak terminal voltage was set to 10,000. vo l t s . The machine was simulated in the steady s tate over f i v e hundred 200-psec. time steps. A l l d-q-0 currents were constant with time. However, when the network was modelled in the Transients Program, very smal l , s tab le o s c i l l a t i o n s about the steady s tate values were observed. The ampl i -As a prel iminary tes t , a very simple network was chosen so that GENERATOR R = 10.fi INFINITE BUS V\AAA VNAAA 0 0 10000.10°\I 5000.A0°V Fig . 9. Simple network conf i gura t ion . 39. t i ldes were very s m a l l , and in the case o f the " d " c u r r e n t , i t amounted to 0.0004 p e r c e n t o f the s t e a d y s t a t e v a l u e . The cause o f t h i s o s c i l l a t i o n c o u l d not be t r a c e d to any e r r o r in the machine program and i t s i n f l u e n c e on the s o l u t i o n was e n t i r e l y n e g l i g i b l e . In the f i r s t e a s e , the s o l u t i o n was a lways a c h i e v e d in one i t e r a t i o n , w h i l e i n t e r f a c e w i t h the T r a n s i e n t s P r o -gram caused conve rgence to be reached a f t e r two i t e r a t i o n s a t each t ime s t e p . 6.2 S i n g l e L i n e - t o - G r o u n d F a u l t The g e n e r a t o r had the same e l e c t r i c a l parameters and was m o d e l l e d as one i n f i n i t e mass as b e f o r e . T h i s t e s t ca se i s t h a t used in Example k o f r e f . [ 5 ] , where a s i n g l e 1 i n e - t o - g r o u n d f a u l t i s a p p l i e d to phase " a " . The network c o n f i g u r a t i o n i s shown in F i g . 10. The s i m u l a t i o n was c a r r i e d out o v e r f o u r hundred t ime s t e p s w i t h a s t e p s i z e o f 500 y s e c . Phase c u r - . r e n t s , f i e l d c u r r e n t and e l e c t r o m a g n e t i c t o r q u e were a lways w i t h i n 0.02 p e r -c e n t o f the r e s u l t s o b t a i n e d by Method I o f [ 5 ] . T h i s s l i g h t d i s a g reement 3 SINGLE PHASE UNITS 13 .8 /130 .8 kV L = 63.2mH 13.95A-30°kV INFINITE BUSBAR ^ TRANSMISSION LINE L = 0.22H L" = 0.096H 1 3 7 - 2 3 / - - 2 0 ° k V FAULTED BUSBAR F i g . 10. Network c o n f i g u r a t i o n f o r s i n g l e l i n e - t o - g r o u n d f a u l t . 40. between the t e s t and the r e f e r e n c e s i m u l a t i o n s was p r o b a b l y caused by a ve ry sma l l d i s c r e p a n c y in the da ta used f o r the two s i m u l a t i o n s . Exac t agreement w i t h Method II o f [5] was a c h i e v e d when i d e n t i c a l t e s t da t a was used in both c a s e s . (The c o m p i l e d v e r s i o n o f Method II was a v a i l a b l e to the au tho r but not Method l>) 6 . 3 IEEE Subsynchronous Resonance Benchmark T e s t Case The g e n e r a t o r i s now a s i x - m a s s r a t h e r than a one-mass mechan i ca l s y s tem, h a v i n g a s e v e n - w i n d i n g e l e c t r i c a l mode l . The IEEE benchmark t e s t ca se f o r the s i m u l a t i o n o f subsynchronous re sonance [22] was based on the a c t u a l e x -p e r i e n c e a t the Mohave g e n e r a t i n g p l a n t in Sou thern Nevada. A s i m p l e RLC c i r c u i t , p r o p e r l y t u n e d , was found to p roduce the same t r a n s i e n t prob lems as those o b s e r v e d in the a n a l y s i s o f the a c t u a l s y s tem. The e l e c t r i c network c o n f i g u r a t i o n i s shown in F i g . 11 and the mechan i ca l model o f the g e n e r a t o r T h i s ca se p r o v i d e s a much more s e v e r e t e s t o f the s o l u t i o n method. B INFINITE BUSBAR X c . - j : v 3 7 1 I— UNITY VOLTAGE c l o s e t=0. open t=0.075s F i g . 11 . E l e c t r i c network f o r the s i m u l a t i o n o f subsynchronous r e sonance . 41 HP ) ' )) |"p \ f)\ LPA ^  J \ L P B \ ^)\QT0R\ ^)) EXC F i g . 12. M e c h a n i c a l model o f the s h a f t s y s tem. i s shown in F i g . 12. The e l e c t r i c a l and mechan i ca l machine d a t a i s g i v e n in r e f . [22]. The machine was s i m u l a t e d in p h y s i c a l u n i t s w h i l e the network was r e p r e s e n t e d in the T r a n s i e n t s Program in the per u n i t s y s tem. T h i s r e -q u i r e d the T h e v e n i n e q u i v a l e n t o f the network from the T r a n s i e n t s Program to be c o n v e r t e d to p h y s i c a l u n i t s a t the b e g i n n i n g o f the t ime s t e p r o u t i n e and then back to per u n i t q u a n t i t i e s a f t e r the s o l u t i o n was f o u n d . A s i m u l t a n e o u s t h r e e - p h a s e f a u l t was a p p l i e d at bus B in F i g . 11 at t ime t = 0 and then removed a f t e r 0.075 s e c . as the c u r r e n t in each phase c r o s s e d z e r o . The p e r i o d o f the s i m u l a t i o n was 0.4 seconds w i t h a t ime s t e p o f 500 y s e c . F i g . 13 shows the s i m u l a t e d c u r r e n t in phase " c " , the s i m u l a t e d e l e c t r o m a g n e t i c t o rque on the g e n e r a t o r r o t o r , and the s h a f t t o rque between r o t o r and e x c i t e r . ( S h a f t t o rque e q u a l s the a n g l e d i f f e r e n c e between the ends o f the s h a f t s e c t i o n m u l t i p l i e d by the s p r i n g c o n s t a n t . ) The growing o s c i l l a t i o n s o f t o r q u e in the l a s t s e c t i o n o f s h a f t would e v e n t u a l l y l e ad to s h a f t f a i l u r e in a r e a l s y s tem. In f a c t , in the a c t u a l c a s e , s h a f t damage d i d o c c u r a t t h i s p o i n t . The r e s u l t s p r e s e n t e d he re a r e in ve ry good a g r e e -ment w i t h those o b t a i n e d by the program d e v e l o p e d by Southern C a l i f o r n i a E d i s o n and P a c i f i c Gas and E l e c t r i c Co. [22], and w i t h the r e s u l t s g i v e n in [5]. V hi. i r 0.16 0.2 0.24 TIME (SEC) F i g . 13. Subsynchronous r e s o n a n c e , phase c c u r r e n t ( u p p e r ) , e l e c t r o m a g n e t i c t o r q u e on the r o t o r ( c e n t e r ) , and E X C - r o t o r s h a f t t o r q u e ( l o w e r ) . 43. 7. EXAMPLES OF PROGRAM APPLICATIONS Examples a re p r e s e n t e d to i l l u s t r a t e some o f the p o s s i b l e a p p l i -c a t i o n s o f the synchronous machine model in the T r a n s i e n t s Program. The same machine i s used in the examples t o p e r m i t a compar i son o f s h a f t t o rques r e s u l t i n g f rom d i f f e r e n t types o f d i s t u r b a n c e s . The g e n e r a t o r and p a r t o f the network from the IEEE benchmark s tudy i s used in t h i s c h a p t e r . 7.1 S u s t a i n e d T h r e e - P h a s e F a u l t f o r an Un loaded G e n e r a t o r The sys tem d iagram shown in F i g . 14 was used to s i m u l a t e a s i m u l -taneous t h r e e - p h a s e f a u l t on an i n i t i a l l y un loaded g e n e r a t o r . A t t ime t = 0, phase " a " t e r m i n a l v o l t a g e was 2 6 . 0 Z 0 0 kV (RMS 1 i n e - t o - 1 i n e ) when a s i m u l -taneous t h r e e - p h a s e f a u l t t o ground o c c u r r e d on the h i gh s i d e o f the D e l t a / Wye t r a n s f o r m e r s . F i g . 15 shows the s i m u l a t e d e l e c t r o m a g n e t i c t o r q u e on the r o t o r and the mechan i ca l t o rque in the s h a f t s e c t i o n between the l a s t low p r e s s u r e t u r b i n e s t a ge (LPB) and the r o t o r . The t ime span o f the s i m -u l a t i o n was 0.5 seconds w i t h a s t e p s i z e o f 300 p s e c . 7.2 F a u l t y S y n c h r o n i z a t i o n o f an Un loaded G e n e r a t o r In t h i s c a s e , a g e n e r a t o r h a v i n g n o - l o a d e x c i t a t i o n and r a t e d t e r m i n a l v o l t a g e i s c o n n e c t e d to the network . For normal s y n c h r o n i z a t i o n , the network and the g e n e r a t o r v o l t a g e phasor s on e i t h e r s i d e o f the c o n -n e c t i o n p o i n t c o r r e s p o n d in magnitude and ang l e at the t ime o f s w i t c h i n g . F a u l t y s y n c h r o n i z a t i o n o c c u r s when the se v o l t a g e s do not c o r r e s p o n d in mag-n i t u d e o r a n g l e . In t h i s s e c t i o n , two f a u l t y s y n c h r o n i z a t i o n s a re s i m u l a t e d . In the f i r s t examp le , the g e n e r a t o r phase " a " v o l t a g e leads the sys tem phase " a " by 120° e l e c t r i c a l ( e l . ) , o r , e r r o r a n g l e 3 = 1 2 0 ° . In the s e c o n d , the e r r o r ang l e i s 9 = 240° e l . In both cases the magni tudes o f the v o l t a g e s kk. 26. Z_0°kV A 4 2 6 . / 5 3 8 . 7 kV SINGLE PHASE UNITS CLOSE AT t=0 CLOSE AT t=0 CLOSE AT t=0 F i g . 14. System d iagram f o r t h r e e - p h a s e f a u l t . o.o o.os o.i o.is 0.2 0.25 0.3 TIME (SEC) 0.35 0.4 0.'45 0.50 F i g . 15. E l e c t r o m a g n e t i c t o r q u e on r o t o r ( u p p e r ) , and L P B - r o t o r mechan i ca l t o rque ( l o w e r ) . h5. at the c o n n e c t i o n p o i n t s were the same. T h i s was a c h i e v e d by a d j u s t i n g the v o l t a g e a t the i n f i n i t e bus so t h a t , w i t h p r o p e r s y n c h r o n i z a t i o n , no c u r r e n t would f l ow between g e n e r a t o r and network. The g e n e r a t o r used p r e v i o u s l y was c o n n e c t e d to the s y s t e m , as shown in F i g . 1 6 . The t r a n s m i s s i o n l i n e to the i n f i n i t e bus was q u i t e l o n g , h a v i n g X Q = 1 . 5 6 pe r u n i t and X ^ = 0 . 5 0 pe r u n i t (based on 8 9 2 . 4 MVA and 5 3 8 . 7 2 kV ) . The f a u l t y s y n c h r o n i z a t i o n s were s i m u l a t e d o v e r a t ime span o f 0 . 5 seconds w i t h a 3 0 0 y s e c . t ime s t e p s . F i g . 17 shows the s i m u l a t e d e l e c t r o -magnet i c t o r q u e and the mechan i ca l t o rque in the L P B - r o t o r s h a f t f o r an e r r o r ang l e o f 120° e l . ( g e n e r a t o r o p e r a t i o n ) . F i g . 18 shows the c o r r e s p o n d i n g to rques f o r an e r r o r ang l e o f 2 4 0 ° e l . (motor o p e r a t i o n ) . These p l o t s com-pare q u i t e w e l l t o those g i v e n in r e f . [ 2 3 ] . Two c o n c l u s i o n s may be drawn from F i g s . 17 and 18. (1) The two f a u l t y s y n c h r o n i z a t i o n s w i t h e r r o r ang l e s o f 1 2 0 ° and 2 4 0 ° e l . , p roduce d i f f e r e n t t o rques - both mechan i ca l and e l e c t r i c a l . (2) A f a u l t y s y n c h r o n i z a t i o n w i t h e r r o r a n g l e 120° p roduces more s e v e r e mechan i ca l and e l e c t r i c a l t o r q u e s than f a u l t y s y n c h r o n i z a t i o n o f 2 4 0 ° e l . These a re e s s e n t i a l l y the same c o n c l u s i o n s reached in [23] f o r a somewhat d i f f e r e n t machine and network . The Brown Bove r i g e n e r a t o r was r a t e d 780 MVA w i t h f i v e masses on the t u r b i n e - g e n e r a t o r s h a f t compared to the IEEE benchmark g e n e r a t o r w i t h an 8 9 2 . 4 MVA r a t i n g and h a v i n g s i x r o t a t i n g masses . 46. 26.^0° kV 3, SINGLE PHASE UNITS 26/538.7 kV TRANSMISSION LINE ZQ=162.6 + J507.3 P. Z ^ . 5 + j 162.6 n 6 6 6 380.95^30 0-6 kV (1) 6=120° (2) 6=240° F i g . 16. System d iag ram f o r f a u l t y s y n c h r o n i z a t i o n . 477' TIME (SEC) Fig. 17. Torques due to faulty synchronization (3 - 120°) TIME (SEC) Fig. 18. Torques due to faulty synchronization (3 - 2k0°) 48. 8. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH T h i s t h e s i s has d e s c r i b e d the e l e c t r i c a l and s h a f t t o r s i o n a l models o f the synchronous mach ine , as i t was added to the U.B.C. E l e c t r o -magnet i c T r a n s i e n t s Program, and a new s o l u t i o n t e c h n i q u e f o r s o l v i n g the comple te s e t o f machine e q u a t i o n s s i m u l t a n e o u s l y w i t h the e l e c t r i c network . The synchronous machine model was i n t e r f a c e d w i t h the T r a n s i e n t s Program by means o f a T h e v e n i n e q u i v a l e n t c i r c u i t o f the network from the machine t e r m i n a l s . The r e l i a b i l i t y o f the s o l u t i o n method has been v e r i f i e d in t e s t s t u d i e s where compar i son was made to the r e s u l t s o f o t h e r e n t i r e l y d i f f e r e n t s o l u t i o n t e c h n i q u e s . Two t e s t s t u d i e s were i n t e n d e d as examples o f p o s s i b l e a p p l i c a t i o n s o f the machine program. In a f a u l t y s y n c h r o n i z a t i o n s tudy i t was p o s s i b l e t o draw c o n c l u s i o n s about the e l e c t r o m a g n e t i c and s h a f t t o r q u e s t h a t agreed w i t h s t u d i e s made by a m a n u f a c t u r e r . The re a re a number o f p o s s i b l e a reas o f a d d i t i o n a l i n v e s t i g a t i o n s . Some o f t he se a r e : (1) i n c l u s i o n o f s a t u r a t i o n e f f e c t s , (2) improvement o f the s t e a d y s t a t e i n i t i a l i z a t i o n o f the e l e c t r i c a l v a r i a b l e s t o i n c l u d e n e g a t i v e and z e r o sequence components in the i n i t i a l i z a t i o n [14], (3) s tudy o f the e f f e c t s on e l e c t r o m a g n e t i c and s h a f t t o r q u e o f mechan-i c a l damping in the s h a f t s y s tem, wh ich i s not y e t b e i n g measured by m a n u f a c t u r e r s , and , (4) e x a m i n a t i o n o f the s h a f t t o rques d u r i n g c o r r e c t and f a u l t y ^ s y n c h r o -n i z a t i o n . The synchronous g e n e r a t o r model s h o u l d p rove to be ve ry u s e f u l in s t u d y i n g s h a f t t o r s i o n a l and o t h e r e f f e c t s in synchronous machines f o l l o w i n g t r a n s i e n t d i s t u r b a n c e s , where the t ime p e r i o d o f s tudy may be as long as one s e c o n d . 50. REFERENCES 1. V. Brandwajn and H.W. Dommel, " I n t e r f a c i n g G e n e r a t o r Models With an E l e c t r o m a g n e t i c T r a n s i e n t s P r o g r a m " , Paper A 76359-0 . .presented a t !the IEEE PES Summer M e e t i n g , P o r t l a n d , O r e . , J u l y 19-23, 1976. 2. IEEE Committee R e p o r t , "A b i b l i o g r a p h y f o r the s tudy o f subsynchronous resonance between r o t a t i n g machines and power s y s t e m s " , IEEE T r a n s . Power App. S y s t . , v o l . PAS-95, pp. 216-218, J a n . / F e b . 1976. 3. W.S. Meyer , T r a n s i e n t s Program Memorandum, December 23, 1977-4. K. C a r l s o n , E.H. L e n f e s t , and J . J . L a F o r e s t , "MANTRAP, Machine and Network T r a n s i e n t s P rog ram" , Paper V -D, p r e s e n t e d a t the 9th PICA C o n f e r e n c e , New O r l e a n s , L a . , June 2-4, 1975. 5. V. B randwajn , " S ynch ronous G e n e r a t o r Models f o r the S i m u l a t i o n o f E l e c -t r o m a g n e t i c T r a n s i e n t s " , Ph.D. t h e s i s , U n i v e r s i t y o f B r i t i s h C o l u m b i a , Vancouver B . C . , A p r i l , 1977-6. W.S. Meyer , T r a n s i e n t s Program Memorandum, March 10, 1974. 7. E.W. K imbark , "Power System S t a b i l i t y : Synchronous M a c h i n e s " , New York : Dover , 1968 ( r e p r i n t o f 1956 e d . ) . 8. B. Adk in s and R.G. H a r l e y , "The Genera l Theory o f A l t e r n a t i n g C u r r e n t Mach ines : A p p l i c a t i o n to P r a c t i c a l P r o b l e m s " , Chapman and H a l l , London, 1975. 9. G. Gross and M.C. H a l l , " Synchronous Machine and T o r s i o n a l Dynamics S i m u l a t i o n in the Computat ion o f E l e c t r o m a g n e t i c T r a n s i e n t s " , to be p u b l i s h e d in the IEEE T r a n s , on Power App. and S y s t . 10. V. Brandwajn and H.W. Dommel, " Synchronous Machine Parameters in A n a l y s i s o f E l e c t r o m a g n e t i c T r a n s i e n t s " , P r o c e e d i n g s o f the Canad ian Communicat ions arid Power C o n f e r e n c e , O c t . 20-22, 1976, pp. 393"395. 11. R.H. P a r k , "Two R e a c t i o n Theory o f Synchronous Mach ines - P a r t I, Gen-e r a l i z e d Methods o f A n a l y s i s " , A l E E T r a n s . , v o l . 48, J u l y 1929, pp. 716-730. 12. Y. Yu , " T h e Torque T e n s o r o f the Genera l M a c h i n e " , T r a n s . Amer. In s t . E l e c t . E n g r s . I l l , v o l . 81, pp. 623-629, Feb. 1963-51. 13- P . J . Fenwick , "Dynamic Response A n a l y s i s o f T u r b o - G e n e r a t o r U n i t s " , U n i v e r s i t y o f B r i t i s h C o l u m b i a , Vancouve r , B .C . , 1976. 14. V. B randwajn , W.S. Meyer , and H.W. Dommel," Synchronous Mach ine \n\t t i a l i z a t i o n f o r Unba lanced Network C o n d i t i o n s W i t h i n an E l e c t r o m a g n e t i c T r a n s i e n t s P rog ram" , to be p r e s e n t e d a t the 1978 IEEE Power E n g i n e e r i n g S o c i e t y Summer M e e t i n g . 15. P.M. Anderson and A .A . Fouad, "Power System C o n t r o l and S t a b i l i t y " , Iowa S t a t e U n i v e r s i t y P r e s s , Ames, Iowa, 1976. 16. H.W. Dommel and N. S a t o , " F a s t 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 T r a n s . Power App. S y s t . , v o l . PAS - 9 1 , J u l y / A u g . 1972, pp. 1 6 4 3 - 1 6 5 0 . 17. G. Gross and A.R. B e r g e n , "A C l a s s o f New M u l t i s t e p I n t e g r a t i o n A l g o -r i thms f o r the Computat ion o f Power System Dynamical R e s p o n s e " , IEEE T r a n s . Power App. S y s t . , v o l . PAS -91, J a n . / F e b . 1977 , pp. 2 9 3 " 3 0 6 . 18. C.W. Gea r , " N u m e r i c a l I n i t i a l Va l ue Problems in O r d i n a r y D i f f e r e n t i a l E q u a t i o n s " , P r e n t i c e - H a l l , Englewood C l i f f s , N. J . , 1974. 19 . H.W. Dommel, " N o t e s on Power System A n a l y s i s " , U n i v e r s i t y o f B r i t i s h C o l u m b i a , Vancouve r , B .C . , 1974. 2 0 . W.S. Dorn and D.D. McCracken , " N u m e r i c a l Methods w i t h F o r t r a n IV Case S t u d i e s " , W i l e y , T o r o n t o , 1972. 2 1 . A . E . F i t z g e r a l d , C. K i n g s l e y , and A. Kusko, " E l e c t r i c M a c h i n e r y " , 3 r d e d . , M c G r a w - H i l l , New Y o r k , 1971. 2 2 . IEEE Subsynchronous Resonance Task F o r c e R e p o r t , " F i r s t Benchmark Model f o r the Computer S i m u l a t i o n o f Subsynchronous R e s o n a n c e " , IEEE T r a n s . Power App. S y s t . , v o l . PAS - 9 6 , S e p t . / O c t . 1977, PP- 1 5 6 5 " 1 5 7 0 . 2 3 . M. Canay, " S t r e s s e s in T u r b o g e n e r a t o r Se t s due to E l e c t r i c a l D i s t u r b -a n c e s " , Brown Bove r i Review, v o l . 6 2 , September 1975, pp. 4 3 5 _ 4 4 3 . 52. APPENDIX 1 TRANSFORMATION OF THEVENIN EQUIVALENT OF THE NETWORK FROM PHASE TO D-Q-0 COORDINATES The Transients Program calculates the three-phase Thevenin equivalent of the network at the synchronous machine terminals, given by equation (19) in phase coordinates. To begin the transformation of (19) to d-q-0 coordinates, equations (3) to (7) must be used. Equation (19) becomes: v 5 d q ° = [ P H Z / ^ H P ] " 1 i s d ^ ° + [P] y with s \$ J l i J s 1 L i J "No (1-1) v. dqo v s d v q V o (1-2) [ Z ^ ^ ] = [VZA;VZB|VZC] ( 1 _ 3 )  VNo = ^ • (1-4) VZA, VZB, VZC and El are column vectors of dimension three. For convenience, they have been given the same alphanumeric names as the corresponding vectors in the Fortran code. The following trigonometric identities w i l l be useful. 2 1 1 cos 0 = -j + Y cos20 (1-5) sin 20 - \ - \ cos20 ( 1 _ 6 ) cos0 sinG = y sin20 ^ 2 2 cos 0 - sin 0 = cos20 (1-8) cos0 cos(0 - 120) = - (1 + cos20) + - | sin20 (1-9) 53. cosG cos(0 + 120) = - f- (1 + cos26) - -™ sin20 (1-10) 4 4 cos0 sin(0 - 120) = (1 + cos20) - 7 sin20 4 4 (1-11) /3 1 cos0 sin(0 + 120) = T (1 + cos2©) - — sin2 0 1 0. 4 4- (1-12) sin0 sin(0 - 120) = ~ (1 - cos20) - 7- sln2e -4 4 (1-13) -1 /3 sin0 sln(©•+'•; 120) .= 7- (1 - cos20) + 7 sin2 0 ,, , .4 4 (1-14) i/3 1 sine cos(0 - 120) = 7- (1 - cos2G) - 7 sin20 1 c. 4 4 (1-15) - / J 1 sin0 cos(0 + 120) = -j (1 - cos20) - j sin20 (1-16) cos 2(0 - 120) =4-7- cos20 - -~ sin20 , 2 4 4 (1-17) 1 i/T cos(0 - 120) sin (0 - 120) = - j sin2 0 + ~ cos2 0 (1-18) 1 /3 cos(0 - 120) sin (0 + 120) = j sin2 0 - (1-19) cos(0 + 120) sin (0 - 120) = y sin2 0 + (1-20) 1 /3 cos(0 + 120) sin (0 + 120) = - f- sin20 - cos20 (1-21) 4 4 cos(0 - 120) cos (0 + 120) = - i+'i--cos2 0 (1-22) 4 2 sin(0 - 120) sin (0 + 120) = - 7 - \ cos20 (1-23) 4 2 cos 2(0 + 120) = j - j cos20 + ~ sin20 (1-24) sin 2 (0 - 120) = j + ~ cos2 0 + sin20 (1-25) sin (0 + 120) = j + i cos2 0 - -| sin20 (1-26) Where 0 is the angle between phase "a" and the direct axis, measured in the direction of rotation. 54. Define [Z 1 by t n e following expression: (1-27) or [ Z P ! 2 3 cosG cos(0 - 120) cos(0 + 120) sinG sin(0 - 120) sin(0 + 120) 1 /2 1 /2 cos 0 1 S2 sin 0 VZA(l) VZA(2) VZA(3) VZB(l) VZB(2) VZB(3) VZC(l) VZC (2) VZC(3) cos(0 - 120) sin(0 - 120) cos(0 + 120) sin(0 + 120) 1 /2 1 /2 1 /2 (1-28) The following can be obtained from (1-28) with equations (1-5) to (1-26) Z (1,1) = -| cos 20 VZA(l) + \ cos0 cos(0 - 120) VZA(2) + -| cos0 cos(0 + 120) VZA(3) + j cos0 cos(0 - 120) VZB(l) + | cos 2(0 - 120) VZB(2) + | cos(0 + 120)cos(0 - 120) VZB(3) + | cos0 cos(0 + 120)VZC(1) + -| cos(0 - 120) cos(0 + 120)VZC(2) + -| cos 2(0 + 120) VZC(3) = (3 + \ c o s 2 0 ) V Z A (D + (~\- \ c o s 2 0 + sin20)VZA(2) 6 + (_ A _ 1 c o s 2 0 _ fl sin20)VZA(3) + (~ \ - \ cos20 + ^ sin20)VZB(l) 55. + ( T - T c o s 2 0 - ^ ? sin29)VZB(2) + (- T + -77 cos20)VZB(3) + \ ~ \ cos29 - i-jr sifL20)VZC(l) + (- -jr + -| cos2G) VZC(2) + ( | - i cos2e + sin29)VZC(3) (1-29) Z (1,2) = ^ sin20 VZA(l) + (- -> sin20 + — - ~ cos20)VZA(2) p 3 6 6 6 + (- i - sin20 - + cos20)VZA(3) + (- |- sin20 - cos20 - ~)VZB(1) + (- f- sin20 + — cos20)VZB(2) + (~„sin20 + ^|)VZB(3) + (- -|- sin20 + cos20 + ~)VZC(1) + ( j sin20 - -^|)VZC(2) 1 V J + (- - sin20 - -|- cos20)VZC(3) (1-30) /2 Z p(l,3) = - j [cos0VZA(l) + cos(0 - 120) VZA(2) + cos(0 + 120) VZA(3) + cos0VZB(l) + cos(0 - 120)VZB(2) + cos(0 + 120)VZB(3) + cos0VZC(l) + cos (0 -120)VZC(2) + cos(0 + 120)VZC(3)] (1-31) Z (2,1) = i sin2@VZA(l) + (- \ sin20 - ^  cos20 - ~ )VZA(2) p J b o b + (- |- sin2© + ^~ cos2 0 + ^)VZA(3) + (- |- sin20 + ^ - ^~ cos20)VZB(l) + (- -|.sin20 + ^ cos20)VZB(2) + (~ sin20 - -5)VZB(3) + (- \ sin20 - + ^ cos2o)VZC(l) o o b o b + (- | sin20 + )VZG/2) + (- | sin20 - cos20)VZC(3) (1-32) 56. Z (2,2) = (4 - T cos20)VZA(l) + (-7 + 7 cos20 - ^ sin20)VZA(2) P J J D O O + ( - 7 + 7 c o s 2 © + sin20)VZA(3) + (- 7 + T> cos20 sin20)VZB(l) 0 0 6 6 6 6 + (4 + T cos20 + ~ sin20)VZB(2) + (- 7 - \ cos20)VZB(3) 0 0 6 6 J 1 1 ' i i + (-7 + 7 cos20 + sin20)VZC(l) + ( - 7 - 7 cos20)VZC(2) 0 0 O O j + ^ + i C O s 2 e ~ "T s i n 2 0 ) v z c ( 3 ) (1-33) /2 Z (2,3) = [sin0 VZA(l) + sin(0 - 120)VZA(2) + sin(0 + 120)VZA(3) P -5 + sin0 VAB(l) + sin(0 - 120) VZB(2) + sin (0 + 120)VZB(3) + sinOVZC(l) + sin(0 - 120)VZG(2) + sin(0 + 120)VZC(3)] (1-34) /2 Z p(3,l) = - j [cos0(VZA(l) + VZA(2) + VZA(3)) + cos(0 - 120)(VZB(1) + VZB(2) + VZB(3)) + cos(0 + 120)(VZC(1) + VZC(2) + VZC(3))] (1-35) /2 Z (3,2) = — [sin (VZA(l) + VZA(2) + VZA(3)) + sin(0 - 120)(VZB(1) + VZB(2) P 3 + VZB(3)) + sin (0 + 120)(VZC(l) + VZC(2) + VZC(3))] (1-36) Z (3,3) = 7 [VZA(l) + VZA(2) + VZA(3) + VZB(l) + VZB(2) + VZB(3) + VZC(l) P 3 + VAC(2) + VZC(3)] (1-37) Define E p by the following expression: E = [P] V p No (1-38) [P ]v N o - [P] E l ( l ) El(2) El(3) (1-39) 57. 3 cosG E l ( l ) + cos(0 - 120)E1(2) + cos(0 + 120)E1(3) sine.El(l) + sln(0 - 120)El(2) + sin(0 + 120)E1(3) 1_ vT /2 /2 El ( l ) + El(2) + El(3) The d-q-0 voltages may now be determined. (1-40) V = Z (1,1) i + Z (1,2)1 + Z (1,3)1 + E (1) Q p d p . q p o p (1-41) or v d = + + + + + + + + + +• + 0j + § cos20)VZA(l) + (- i - i - cos20 + 3 sin20)VZA(2) - i - -|- cos20 - ~ sin20)VZA(3) - T - T cos20 + sin20)VZB(l) O D D ~ - j cos20 - ~ sin20)VZB(2) + (-\ + \ cos20)VZB(3) ~ \ ~ \ cos20 - ~ sin20)VZC(l) + (- |- + -j cos20)VZC(2) j - |- cos20 + sin20)VZC(3)}- i d sin20VZA(l) + (- |- sin20 + - cos20)VZA(2) - I - sin20 - Sr + cos20)VZA(3) 6 6 6 1 /3 / J sin20 - -~ cos2G - — )VZB(1) 6 6 - \ sin20 + ^ | cos20)VZB(2) + (^  sin20 + )VZB(3) b o 5 o _ 1 s i n 2 9 + .i f c o s 2 0 +3. )VZC(1) + (| sin20 - 3)vzc(2) 6 6 6 1 /3 - sin20 - ~ cos20)VZC(3) }• i 6 6 q /I + ~ {cos0 VZA(l) + cos(0 - 120)VZA(2) + cos(0 + 120)VZA(3) 58. + cos© yZB(l) + cos (.9 - 120)VZB(2) + cos,(.0 + 120)VZB(3) + cose VZC(l) + cos(0 - 120)VZC(2) + cos(0 + 120)VZC(3)}•! [cos0E 1(1) + cos(0 - 120)E1(2) + cos(0 + 120)E1(3)J 0 (1-42) Equation (1-42) can be rewritten as: • v. = cos20[2VZA(l) - VZA(2) - VZA(3) - VZB(l) - VZB(2) + 2VZB(3) d o - VZC(l) + 2VZ(2) - VZC(3)] + ~ sin20[VZA(2) - VZA(3) + VZB(l) - VZB(2) - VZC(l) + VAC(3)] + \ [2VZA(1) - VZA(2) - VZA(3) - VZB(l) + 2VZB(2) - VZB(3) o - VZC(l) - VZC(2) + 2VZC(3)]} i , d + { i sin20[2VZA(l) - VZA(2) - VZA(3) - VZB(l) - VZB(2) + 2VZB(3) b - VZC(l) + 2VZC(2) - VZC(3)] /3 + ~ cos20[-VZA(2) + VZA(3) - VZB(l) + VZB(2) + VZC(l) - VZC(3)] o + - f [VZA(2) - VZA(3) - VZB(l) + VZB(3) + VZC(l) - VZC(2)]} i b q + - j {cos0[VZA(l) + VZB(l) + VZC(l)] + cos(0 - 120) [VZA(2) + VZB(2) + VZC(2)] + cos(0 + 120)[VZA(3) + VZB(3) + VZC(3)]} i + [ c o s 0 E l ( l ) + cos(0 - 120)E1(2) + cos (0 + 120)El(3) ] (1-43) The same procedure for V and V gives: q o vq • y ^ • h + y 2 > 2 ) • \ + z p ( 2 > 3 ) • h + E P ( 2 ) of, (1-44) 59. = {~ sin2Gl2VZA(l) - VZA(2) - VZA(3) - VZB(l) - VZB(2) + 2VZB(3) - VZC(l) + 2VZC(2) - VZC(3)] /3 + -g- cos20[-VZA(2) + VZA(3) - VZB(l) + VZB(2) + VZC( l ) - VZC(3)] /3 - - j [VZA(2) - VZA(3) - VZB(l) + VZB(3) + VZC( l ) - VZC(2)]} i + {- j COS20I2VZA(1) - VZA(2) - VZA(3) - VZB(l) - VZB(2) + 2VZB(3) - VZC ( l ) + 2VZC(2) - VZC(3)] + |- [2VZA(1) - VZA(2) - VZA(3) - VZB(l) + 2VZB(2) - VZB(3) - VZC(l) - VZC(2) + 2VZC(3)]} i q /2 + ~ {siri0[VAZ(l) + VZB(l) + V Z C ( l ) ] + sin(0 - 120)[VZA(2) + VZB(2) + VZC(2) + sin(0 + 120)[VZA(3) + VZB(3) + VZC(3)]} i vo = V3'1* ' + Z P ( 3 ' 2 ) * \ + Z P ( 3 ' 3 ) * "0 + E p ( 3 ) /2 VQ = {cos0[VZA(l) + VZA(2) + VZA(3)] + cos(0 - 120)[VZB(1) + VZB(2) + VZB(3)] + cos(0 + 120)[VZC(1) + VZC(2) + VZC(3)]} i £ + — {sin0[VZA(l) + VZA(2) + VZA(3) ] '-+ sin(0 - 120)[VZB(1) + VZB(2) + VZB(3)] + sin(0 + 120)[VZC(1) + VZC(2) + VZC(3)]} i (1-46) 60. + | (VZA(l) + VZA(2) + VZA(3) + VZB(l) + VZB(2) + VZB(3) + VZC(l) + VZC(2) + VZC(3)} i + - [El(l) + El(2) + El(3)] ^ (1-47) The following alphanumeric names are used in the code and are defined by: CONSAl ==-— (VZA(l) + VZB(l) + VZC(l)) CONSA2 = ^  (VZA(2) + VZB(2) + VZC(2)) C0NSA3 = ^ | (VZA(3) + VZB(3) + VZC(3)) CONSA4 = VZA(l) + VZA(2) + VZA(3) C0NSA5 = VZB(l) + VZB(2) + VZB(3) CONSA6 = VZC(l) + VZC(2) + VZC(3) CONI0C = -| (C0NSA4 + C0NSA5 + C0NSA6) AC0N1A = 2VZA(1) - VZA(2) - VZA(3) - VZB(l) - VZB(2) + 2VZB(3) - VZC(l) + 2VZC(2) - VZC(3) AC0N1B = VZA(2) - VZA(3) + VZB(l) - VZB(2) - VZC(l) + VZC(3) A 1 / J ACONSl = y cos20 • AC0N1A + — r sin20 • ACONIB 6 6 ACONS2 = T> sin20 • ACONIA - — cos20 • ACONIB 6 6 BCONS1 = \ ' ACONIA + ~ (VZB(2) - VZB(3) - VZC(2) - VZC(3)) D Z BCONS2 = — (VZA(2) - VZA(3) - VZB(l) + VZB(3) + VZC(l) - VZC(2)) CONI0A = cos0 • CONSAl + cos(0 - 120) • CONSA2 + cos(0 + 120) • CONSA3 CONSTA = C0NI0A • i + ~ (cos0El(l) + cos(0 - 120)E1(2) o 3 + cos(0 + 120)E1(3)) 61 . A CONIOB = sin0 • CONSA1 + sin(0 - 120) • C0NSA2 + sin(0 + 120) • CONSA3 CONSTB = CONIOB • i + t./- (sin0 • E l ( l ) + sin(0 - 120)E1(2) o V J + sin(0 + 120)E1(3)) DCONS1 = (cos0 • CONSA4 + cos(0 - 120) • CONSA5 + cos(0 + 120) • CONSA6) DC0NS2 = — (sin0 • CONSA4 + sin(9 - 120) • CONSA5 + sin(0 + 120) • CONSA6) Taking advantage of the symmetry in (1-43) and (1-45), the expressions for V,, V and V reduce to: d q o V d = (AC0NS1 + BC0NS1) • i d + (AC0NS2 + BCONS2) • i + CONSTA (1-48) V = (AC0NS2 - BCONS2) • i , + (-AC0NS1 + BC0NS1) • i + CONSTB q d q (1-49) V = DC0NS1 • i , + DCONS2 • i + CONIAC • i 0 d q o + - (El(l) + El(2) + El(3)) /3 J (1-50) The partial derivatives of v d » V q » ^ 0 with respect to the electric angle 9 may be found quite readily from (1-41), (1-42) and (1-43). 3V • ~ = -2 • AC0NS2 • i , + 2 • AC0NS1 • I - CONSTB 30 d q v . l -31; 3 V ^ - = 2 • AC0NS1 • i , + 2 • AC0NS2 • i + CONSTA (1-52) ofc) d q 3Vn -rrp- = - DCONS2 • i , + DC0NS1 • i (1-53) 30 d q 62. From (1-48), (1-49), and (1-50), the p a r t i a l d e r i v a t i v e s w i t h r e s p e c t t o i d , i q , i Q a r e : T - r 2 - = AC0NS1 + BC0NS1 0"54) d -zA = AC0NS2 + BC0NS2 " (1-55) o I q -^A = C0NI0A (1-56) 9 , o 8v ^rA = AC0NS2 - BC0NS2 (1-57) 9 l d 9v ^A = - AC0NS1 + BC0NS1 (1-58) o I q av •TA = CONIOB (1-59) a'o — ° = DC0NS1 (1-60) 8 , d av ^ = DC0NS2 (1-61) - r ^ a i q 3v 0 = C0NI0C (1-62) 9 ,o APPENDIX I I STEADY STATE PROGRAM "MSTEA" I TERMINAL SYSTEM FORTRAN G ( 4 1 3 3 6 ) MSTEA 64. SUBROUTINE M S T E A ( I , OMEGA, DELTAT ) COMMON / V L A O / C M ( 7 , 7 ) , C T ( 3 , 3 ) , C I F ( 1 0 ) , Z T R ( 7 0 ) t Z T X { 7 0 ) , Z R < ? . 0 ) CCMMON /VL AD/ V Z A ( 3 0 ) , V Z R ( 3 0 ) , V Z C ( 30 ) ,VZD ( 3 0 ) , E H 3 ) , C U R R T { 3 ) COMMON / V L A O / THET<10J ,ZX<20 I , Z M L ( 9 0 ) , Z M R ( 5 C ) , T E M ( 1 0 ) , A N G , Z 3 COMMON /VL AD/ M H 1 0 ) , M A C H COMMON /MACHIN/ C C N S T 1 , C O N S T ? , C C N S T 3 , H C T , H D T I , A C O N COMMON /MACHIN/ AC0N1,AC0N2,AC0N3, A C O N 1 4 , Z T Y { 7 0 ) , L P A R K { 7 , 7 ) COMMON /MACHIN/ R INEPT181 ,ST I F F ( 7 ) ,DAMp (8) , TORQUE (7) ,THEOLD(8 ) COMMON /MACHIN/ WOLD I8),KBAND(7,7),HIST1(8),HIST2(8 I COMMON /MACHIN/ NPOLES,NMASS,NROTOR DOUBLE PRECIS ION C M , C T , C I F , Z T R , Z T X , Z R , V Z A , V Z B , V Z C , V Z D , E 1 , C U P R T , 1 T H E T , Z X , Z ^ L f Z M R , T E M , A N G . Z 3 DOUBLE PRECIS ION CONST 1, C0NST2, CON ST 3, HOT , H-CT I , ACO N,ACON1.AC0N2, 3 A C O N 3 , A C 0 N 1 4 , Z T Y , L P A R K , R I N E R T . S T I F F , D A M P , T O R Q U E , T H E O L D , W O L D , 3KBAND,H I ST1 ,H I ST2 REAL *8 A S , B S , C S , 0 S , F S , A R , B R , C F , A I , B I , C I , A , B , C , D , E , F , A 1 , B 1 , 3 AL, B L , C l ,DL , D E L T A T , OMEGA,T E, AC0K8, AC0N9,ACCN10.ACONI 1.ACON12, SAC0N13,AC0N1A,ACONIB ,A CONS 1 ,ACON2A,ACON2R,ACCNS2,BCDNS1,PCONST, 3THMECH,WMECH,G(8 I ,DT (7 ,7 I ,RZ (8 I . D E T . P P C L E S t T I M E . T O R Q M INTEGER IPERM(14) ,PARA Ml,PARAM2,PARAM3 ISAVE^I N = U *2> /3 WR ITE (6 ,600 ) N 600 FORMAT ( IH , 20X, • SOLUT ION FOR MACHINE NUMBER ' , 1 4 ) MI(N)=7 I F ( Z M L ( N * 9 - 2 ) . EQ .O .DO ) M H N ) = 6 MU1 = MI { N ) PARAMl=NMASS- l PARAM2 = NR0T0R- 1 PCONST=2.DO/DFLOAT(NPOLES) PP0LES=NP0LES/2 I F ( N . G T . l ) GO TO 610 HDT=0ELTAT/2 .D0 HOT I=2 .DO/DELTAT A N G = ( ? . D O / 3 . D O ) * 3 . 1 4 1 5 926535^00 C G N S T 1 = D S C R T ( 3 . 0 C » CCNST2=CCNST1/DSQRT< 2 .00) CCNST3=0.500*CONST1 ACON=DS0RT (2 .O0 ) /3 .D0 AC0NL=1.D0/CONST2 AC0N2=1.D0/6 .D0 AC0N3=C0NST1*AC0N2 AC0N14=1.DO/DSQRT(2.DO) 610 ILM^N*5-4 ILN=N*2-1 ZMR{ ILM) = ZMR{ ILM|4ZR( ILN) + 2.D0*ZR( ILN+1) ILM=9*N AS=ZML<ILM)*OMEGA CS=ZVL( ILM-8)*CMEGA DS=ZML(ILM-7) FS=ZML( ILM-4)*0MEGA ! L M = 5 * ( N - l ) BS = ZMR( ILM + 1) C CALCULATE POS IT IVE SEQUENCE CURRENTS * * * * * * * * * * * * * * * * * * * * * * * * J = I + 1 K=I + 2 AR=VZA(I) BR=VZA(J) TERMINAL SYSTEM FORTRAN G ( 4 1 3 3 6 ) MSTEA 65. CR=VZA(K) AI = V Z B(II BI= V Z B (J) CI=VZR(K) A M 3 R + C R ) * 0 . 5 D 0 B=(C I -B!>*CONST3 C=(C I4-B I ) *0 .5D0 D=(BR-CR)*CONST3 BR=(AR*BR*CR)/3.DO BI={AI+RI+CI ) /3.DO CR=( AR-,A-eJ/3.D0 C I = ( A T - C - n ) / 3 . D O AR={AR-A+B)/3.D0 AI=(A I - C * D 1 / 3 . 0 0 W R I T E ( 6 , 6 1 l » A R » A I 611 FORMAT{• « , " P O S I T I V E SEQUENCE CURRENT: REAL PART= « , D 1 5 . 7 , 5 X , S ' IMAGINARY PART= ' , 0 1 5 . 7 ) WR ITE (6 ,612 ) 612 FORM AT(1H • ' INFORMATION ABOUT THE UNBALANCE OF THE SYSTEM «) WR ITE (6 ,613 ) 8R,B I 613 FORM A T ( ' » , » Z E R O SEQUENCE CURRENT: REAL PA RT = ' . D 1 5 . 7 . 5 X , a»IMAGINARY PART= « , D 1 5 . 7 > WR ITE(6 ,61A) CR,C I 614 FORMAT{• « , ' N E G A T I V E SEQUENCE CURRENT: REAL PART=' » , D 1 5 . 7 , 5 X , 3 ' IMAGINARY PART = % D 1 5 . 7 > : CALCULATE POS IT IVE SEQUENCE CURRENTS IB ANC IC * * * * * * * * * * * * * * * * BR=-0.5D0*AR+C0NST3*A I BI= - C . 5 D 0 * A I - C C N S T 3 * A R CR=-0.5 00*AR -CCNST3*4 I CI= - 0 . 5 D C * A I * C C N S T 3 * A R C FIND THE ELECTR ICAL POSITION OF THE ROTOR * * * * * * * * * * * * * * * * * * * * MN=N*2 C L = Z R I M N - 1 ) - Z R ( M N » 0L=ZX (MN-1 ) -ZX (MN) AL=ZR(M\> * ( A R+BR + CR) BL=ZX(MN)*( AI + BI+CI ) A L - - A L + B L VZC ( I ) = V Z C { I » - A L - C L * A R * D L * A I V Z C < J ) = V Z C ( J ) + A L - C L * B R * D L * B I VZCC K)=VZC(K) + A L - C L * C R + D L * C I AL=ZR(MNI * {A I *B I *C I I BL=ZX(MN)*(AR+BR+CR) AL=AL+BL vzn{ n = vzo( D - A L - C L * A I -OL*AR VZD( J)=VZD( J ) - A L - C L * B I - D L * B R VZDt K ) = V Z 0 1 K ) - A L - C L * C I - D L * C R A=VZC( I ) B=VZDI I ) C=A*BS*AR-FS*A I 0=B+BS*A I*FS*AR . C=DATAN2(D,C l C ELECTR IC ANGLE OF THE R O T O R * * * * * * * * * * THET(N )=C+1. 5707<5632679500 C FIND I D , I C » V D , V Q TC CALCULATE THE CURRENT IF * * * * * * * * * * * * * * * * * * * HM=N*7-6 A1=DSQRT(A**2+B**2 ) B i=DATAN2(B» A) B=B1-C TERMINAL SYSTEM FORTRAN G(41336> MSTEA 66. F = C O N S T 2 * A l * D S I N ( B ) ZTX l I LM )=F B=CONST2*A l *OCOS(B ) ZTX{ILM+1)=B A1=DSQRT(AR**2+A1**2 ) B l = O A T A N 2 ( A ! « A R ) A=B1-C D=CONST2*Al*DS INC A) ZTR( ILM)=D A=CCNST2*A l *OCOS (A ) ZTR( ILM+1)=A ILM=ILM+2 ZTR ( I LM)=0 .00 ZTX( ILM )=O.DO K=MI(N)-3 DO 620 L=1,K J=L+ILM Z T R U I = 0 . 0 0 620 Z T X ( J ) = 0 . D 0 C= (B+8S*A -CS*D ) / (DS*0MEGA? Z T X ( I L M ^ 1 I = - Z M R ( ( N - l > * 5 * 2 ) * C K=ILM+1 ZTR(K)=C C FORM INDUCTANCE M A T R I X * * * * * * * * * * CALL I N 0 U ( Z M L » C 0 N S T 2 , L P A R K ,N,MU1) I LM=(N-1 ) *7 DO 621 J = 1,M|J1 ILN=ILM+J Z T Y ( T L N I = 0 . 0 0 DO 622 K=1,MU1 62? ZTY( ILN)=ZTY( ILN )+LPARK ( J , K J *ZTR (K > 621 CONTINUE C TRANSFORM VOLTAGES AND CURRENTS BACK TO COIL VARI A B L E S * * * * * * * * * * TE=TMET(N I AC0N8=DC0SfTE) AC0N9=DCCS(TE-ANG) AC0N10=DC0S(TE-ANG) AC0N11=0S IN (TF) ACCN12=DS IN(TE-ANG) AC0N13=DSIN(TE+ANG) AC0NIA=AC0N1*(AC0N8*D+AC0N11*A) AC0NIB = ACCN1* (ACC fv 9* 0* ACON12*A I AC0NS1=ACGN1*(ACON10*D+ACON13*A) AC0N2A=ACCN1*{ACCN8*F+AC0N11*BJ AC0N2B=AC0N l * (ACCN9*F*ACON12*B ) AC0NS2=ACCNl* ( AC0N10*F +AC0N13*B I • TE = P P 0 L E S * ( Z T Y U L M + 1 >*A-ZTY( I L M * 2 » * D ) W R I T E < 6 , 6 2 3 » 0 , A , C , F , B , T H E T < N J 623 FORM AT I • « , ' I - D = « ,D1 5 .7 , 8 X, • I-C= • , 01 5. 7 , E X, • I- F= « , D 1 5 . 7 , / , S « V - D = » , 0 1 5 . 7 , 8 X , « V - Q = • , D 1 5 . 7 , 8 X , « E L E C ANGLE (RAD IANS)= ' , 0 1 5 . 7 ) C I N I T I A L I Z E THE MECHANICAL EQUATIONS * * * * * * * * * * THMECH=PCCNST*THET{N) VMECH=PCCNST*OMEGA WR ITE (6 ,624J TE ' 6 2 4 FORMAT(* • , « E L E C T R C M E C H A N I C A L TCRQUE: R O T O R = • , D l 5 . 7 , 5 X ) I F ( N M A S S . G T . l ) GO TO 650 B=0.D0 RTNERTt 1 ) = U D 1 0 4 TERMINAL SYSTEM FORTRAN G141336) MSTEA 67. OAMP{11*0.00 TnPOUF( l)=TE THEOLO(1)=THMECH WOLDl1)=WMECH G( 1 ) =THMECH GO TO 682 650 I LM= (N -1 ) *7 TORO M=T E IF INROTOR.EQ.NMASS ) TO ROUE(NROTOR 1 = 0 . 0 0 I F fTO&OUEtNPOTOR) .EQ .O .DO) GO TO 630 TORQUE INROTOR>=-( PPOLES*C ) * {ZTXl ILM+4)+TORQUE <NROTOR)*C)/WMECH T0RQM=TE-T0R0UEINPCT0R) WRITE16,625) TORQUE(NROTOR) 625 F O R M A T l • + « , « E X C I T E R = ' , D 1 5 . 7) 630 DO 632 1=1,NM4SS 632 TORQM=TORQM+WMECH*CAMP(I} DO 635 I=1,PARAM2 635 TORQUE( I )= (TORQUE(11/100 .D0 ) *TORQM WRIT E( 6 , 6 2 6 ) ( 1, TORQUE ( I ) , I=1,PARAM2) 626 FORMAT! ' • , • MECHANICAL TORQUES FROM T U R B I N E S * » / » 7 I * ( • , 11»* ) ' » 3 D 1 3 . 6 , 2 X ) ) IF (NROTOR.EQ.NMASS) GO TO 652 A=DA MP(NROTOR) B=RINERT(NROTOR) ILN=NR0TCR*1 DO 651 I=ILN,NMASS J=I -1 OAMPlJ ) = CAMP( I | 651 R I N E p T ( J ) = R I N F P T { I ) CAMPt NMASS) = A RINERTINMASS)=B 652 DO 655 I=1,PARAM1 DO 655 J = l , P A R A M l 655 KBAND l I , J ) =O .DO DO 656 I=1,PARAM1 656 KBANDII , H = S T I F F ( I > I F f N R 0 T 0 R . L E . 2 1 GO TO 658 CO 657 I=2,PARAM2 J = I- 1 K B A N D l I , I )=KBA N D ( I , I ) + ST I F F { J ) K B A N D ( J , I ) = - S T I F F I J ) 657 KBANDlI , J ) = KBANOtJ , I> IF INR0T0R + 1.GE.NMASS ) GO TO 660 658 I F INMASS .EQ.2 ) GC TO 660 ILN=NMASS-2 DO 659 I = NROTOR, ILN J=I*1 K B A N D l I , I ) = K B A N D l 1 , 1 ) + S T I F F ( J ) KBANDl J , I ) = - S T I F F U ) 659 K B A N D d , J ) = KEAND IJ , I) 660 DO 661 I=1,PARAM1 661 G ( I I =TOPQUE l I ) -DAMP( I ) *WMECH DO 662 1=1,NMASS 662 WOLD(I ) = WMECH IF INROTOR.EQ.NMASS I GO TO 664 GlNROTOR)=GINROTOR)+STIFF INROTOR)*THMECH I F INROTOR.EQ.1 ) GC TO 665 664 I=NR0T0R-1 M TERMINAL SYSTEM FORTRAN GH1336) MSTEA 68. 6 ( I ) a R ( I J * S T l F F ( I)*THWECH 665 DET=1.D-14 CALL DSLIMPf KBAND.OT, G, THEOL D, R Z , I PER M , PAR AM 1, 7, DE T, 1, 14) WRITE(6,663) (TH ECLD(I ),1 = 1,PAP AMI ) 663 FORM AT(• * , ' T H E T A - S / S ' ,5X,7(D15.7) ) THEOLD(NMASS)=THMECH IFtMROTOR.FO.l) GO TO 670 DO 669 I=1,PARAM2 669 G(I)=THECLD(I) 67C G(NROTOR)=THMECH IF<NROTCR.EC.NMASS) GO TO 672 ILN=NR0TCR-1 DO 671 I=ILN,NMASS 671 G( I )=THEOLD(I - l ) C BEGIN CALCULATION OF HISTORY TERMS FOR MECHANICAL PART * * * * * 672 B=0.D3 C=O.DO D=G11) DO 675 I=1,PARAM1 E = STIFF<n F=G( I + l ) HIST21I )=-B*C+(B+E)*D-E*F B=E C = D 675 D=F HIST21NMASS)=-B*C*B*D B=HIST2{NROTOR ) IFtNROTOR.ECNMASS) GO TO 677 DO 676 J=NPOTOP,PARAM1 676 H I ST2U J=HIST2I J + l ) 677 UM=(N-1)*7 682 HIST2lNMASS)=B+PP0LES*tZTY(ILM+l)*ZTR{ILM+2)-ZTY(ILM+2)*ZTR(ILM+1 3) 00 678 J=1,NP0TCP 678 H IST2U)=HIST2(J) -TORQUF(J) IF< NMASS.GT.NROTCR ) HI ST2{NROTOR) = HI ST 2(NROTOR)-TOROUE(NROTOR) DO 679 I=1,NMASS HIST2(I)=HIST2(I)*HOT/RINERT{I) 679 HIST1{I)=(HDTI+DAMP(I)/RINERT(I))*THEOLD{I)+2.D0*WMECH 68C ILM=(N-1I*7 THET(NI*THET(N)-1.5707963300 A=O.DO B=O.DO C=O.DO IFCNMASS.EQ.l> GC TO 681 B=(GtPARAM2)-G(NROTOR) )*ST IFF<PARAM2) IF(NMASS.EQ.2) GC TO 681 A=(G(NROTOR-2)-G(PARAM2))*STIFF tNROTOR-2) IFtNMASS.GT.NROTOR) C=(G(NROTOR)-G(NROTOR+l> )*STIFF{NROTOR) 681 WRITE18) AC0N1 A»ACCN1B,AC0NS1» ACON 2 A, AC0N2B»AC0NS2, 3ZTR( ILM+4),TE,A,B,C,THET(N) WRITE(6,673) NVASS,N 673 FORMAT!• »,» INITIAL MECHANICAL ANGLES OF MASSES 1 TO ' , 1 1 , • FOR " 3CHINE ' , 1 2 ) WRITE(6,674) { I,G( I),I=1,NMASS) 674 FORM A T ( ' ' ,4(•MASS-• ,11,D I5.7,5X)1 I=ISAVE RETURN APPENDIX I I I TIME STEP PROGRAM 1 M T E R M I N A L S Y S T E M F O R T R A N 0 ( 4 1 3 3 6 } GEN 70. S U B R O U T I N E G E N U , T I M E , D E L T A T , C M E G A ) COMMON / V L A O / CM ( 7, 7 ), CT{ 3 ,3 ),C IF ( 10 ) , Z TR ( 7C ) , Z T X ( 70 ) , ZR ( 2 0) COMMON / V L A O / V Z A ( 3 0 ) , V Z B ( 3 0 ) t V Z C I 3 0 l , V Z D ( 3 0 ) , E l ( 3 ) , C U R R T ( 3) COMMON / V L A O / T H E T ( 1 0 ) , Z X ( 2 0 ) , Z M L ( 9 0 ) , Z M R < 5 0 ) , T E M < 1 0 ) , A N G , Z 3 CCMMCN / V L A D / M I ( 1 0 ) , M A C H COMMON / N"ACH IN / C O N S T l , C 0 N S T 2 , C C N S T 3 , H D T , H D T I , A C ON COMMON / M A C H I N / ACON 1 , A C O N 2 , A C O N 3 , A C O N 1 4 , Z T Y ( 7 0 ) , L P A R K ( 7 , 7 ) COMMON / M A C H IN/ R I N E R T ( 8 ) » ST I F F ( 7 ) » D A M P ( 8 ) » T O R Q U E { 7 1 t T H E O L D ( 8 ) COMMON / M A C H I N / W O L D ( 8 ) , K B A N D ( 7 , 7 ) , H I S T 1 I 8 ) , H I S T 2 ( 8) COMMON / M A C H I N / N P O L E S , N M A S S , N R O T O R D O U B L E P R E C I S I O N C M , C T , C I F , Z T P , Z T X , Z R , V Z A , V Z B , V Z C , V Z D , E l . C U R R T , 1 T H E T , 7 X , Z M L , Z M R , T E M , A N G , Z 3 D O U B L E P R E C I S I C N CONST I , C O N S T 2 , C O N ST 3 , H D T , H O T I , A C G N , A C C N 1 , A C ON 2, 3 A C 0 N 3 , A C C N 1 4 , Z T Y , L P A R K , R I N E R T , S T I F F , D A M P , T C R Q U E , T H E O L D , W O L D , a K B A N D , H I S T l , H I S T 2 REAL * 8 W , T H E , T H F 2 , D E L T A T , 0 M E G A , A C 0 N 4 , A C 0 N 5 , A C G N 6 , A C 0 N 7 , 3 X 0 L P I 7 ) , V 0 L D ( 7 ) , P S I 0 L 0 ( 7 ) , X ( 7 ) , V ( 7 ) , P S I 1 7 ) , C ( 1 5 ) , C C N S A 1 , C O N S A 2 , 3 C 0 N S A 3 , C O N S A 4 , C O N S A 5 , C O N S A 6 , C 0 N I O C , R V E C ( 7 ) , J A C O B ( 1 5 , 1 5 ) , B C 0NS2 , SCON I O A , C O N I O B , D C C N S 1 , D C O N S 2 . C 0 N S T A , C O N S T B , C A ( 1 5 , 1 5 ) , D E L T A X I 1 5 ) , 3 R Z U 5 ) , D E T . T E , A C 0 N 8 , A C 0 N 9 , AC ONI 0 , ACON 1 i , A C C N 1 2 , A C 0 N 1 3 , A C O N I A , 3 A C 0 N 1 B , A C C N S 1 , A C O N 2 A , A C 0 N 2 B , A C O N S 2 , B C O N S 1 , T H E N E W ( R ) , T H E T A S ( 8 ) , 3 P P 0 L E S , K N E W Q T , B , C , D , E, F , T I M E , HN E WI ( 8 ) , H N E W 2 ( 8 ) I N T E G E R I P E R M { 3 0 ) , P A R A M 1 , P A R A M 2 , P A R A M 3 N = ( I + 2 ) / 3 I S A V E = T MU1=MI (N) P A R A M U K M A S S - 1 P A R A M ? = N P 0 T 0 R - 1 PARAM3=NMASS+MU1 P P 0 L E S = N P C L E S / 2 I T E R = 0 ICONV=0 DO 5 9 3 1 = 1 , N M A S S 5 9 3 THENEW{ I ) = T H E O L D ( I ) + W O L D ( I ) * D E L T A T ! L M = { N - 1 ) * 7 0 0 5 9 4 1 = 1 , M L 1 X O L D U ) = Z T R ( ILM+I ) X( I ) = XOL DI I ) V O L D ( I ) = Z T X f I L M + I } 5 9 4 P S I O L D ( I ) = Z T Y ( I L M + I ) DC 5 9 5 1 = 1 , 5 I L M = N * 5 - I + 1 I L N = 7 - I + 1 595 RVEC U L N ) = Z M R { I L V ) I L N = 2 * N R V E C ( 2 ) = R V E C ( 3 ) - 3 . D 0 * Z R < I L N ) R V E C ( i ) = R V E C ( 2 ) I 1 M S A V E I 2 = I S A V E + 1 13 = 1 S A V E + 2 , C 0 N S A 1 = - A C 0 N * { V Z A ( I 1 ) + V Z B ( I 1 ) + V Z C ( I D ) C C N S A 2 = - A C 0 N * ( V Z M I ? ) * V Z B ( 12 ) * V Z C f 12 I t C O N S A 3 = - A C O N * ( V Z A U 3 ) + VZB ( I 3 ) + V Z C ( 13 ) ) CCNS A 4 = - ( V Z A ( I 1 ) + V Z A ( I 2 ) +V ZA(13 ) ) CONS A 5 = - ( V Z B (I I ) + V Z B ( I 2 ) *V ZB 113 I > C 0 N S A 6 = - ( V Z C ( I 1) +VZCI I 2 ) + V Z C U 3 ) ) C O N I 0 C = ( 1 . D 0 / 3 . 0 0 ) * ( C 0 N S A 4 + C 0 N S A 5 + C 0 N S A 6 ) A C 0 N 1 A = - 2 . D 0 * V Z A U 1 ) + V Z A ( I 2 ) + V Z A ( I 3 ) + V Z B ( 1 1 ) + V Z B ( 1 2 ) - 2 . D 0 * V Z B ( 13 71. ^ T E R M I N A L S Y S T E M F O R T R A N G 1 4 1 3 3 6 ) GEN 3 + V Z C ( T 1 ) - 2 . D 0 * V Z C ( 1 2 >*VZC( 13) AC0N1B = - V Z A ( 1 2 ) + V Z A ( 1 3 ) - V Z B ( I 1 ) + V Z E ( I 2 ) * V Z C < I 1 ) - V Z C ( 1 3 ) I F f N M A S S . E O . l I GO TO 6 0 1 • C * * * B E G I N C O N S T R U C T I O N OF J A C O B I AN M A T R I X - M E C H A N I C A L P A R T F I R S T * * * C * * * U P P E R L E F T S O U A R F ( T I M E I N D E P E N D E N T ) * * * DO 5 9 6 J = 1 , P A R A M 1 DO 5 9 7 K = 1 1 P AR AM I 5 9 7 J A C O B ( J , K ) = K B A N C ( J , K ) * H O T / R I N E R T ( J ) 5 9 6 J A C O B ( J , J ) = J A C 0 3 ( J , J ) + H D T I + D A M P ( J ) / R I N E R T ! J ) C * * * U P P E R R I G H T P A R T ( T I M E I N D E P E N D E N T ) * * * DO 5 9 3 J = l , P A R A N i DO 5 9 8 K= N M A S S , P A R A M 3 5 9 8 J A C O B ( J t K ) = 0 . D O C = O . D O B = S T I F P ( P A R A M 2 ) J A C O B ( P A R A M 2 , P A R A M 3 ) = - H O T * B / R I N E R T ( P A R A M 2 ) I F ( N R O T O R . E O . N M A S S ) GO T O 5 9 9 C=ST I F F ( N R O T O R ) J A C O B ( N R C T 0 R , P A R A M 3 ) = - H D T * C / R I N E R T ( N R O T O R ) C * * & L O W E R L E F T P A R T , S T I L L T I M E I N D E P E N D E N T - E X P L O I T S S Y M M E T R Y * * * 5 9 9 DO 6 0 2 I = N M A S S , P A R AM 3 0 0 6 0 2 J = 1 , P A R A U 1 6 0 2 J A C O ^ U , J ) = J A C O B ( J , I ) C * * * D I A G O N A L E L E M E N T ON LOWER R I G H T CORNER OF J A C O B I A N * * * 601 J A C O B l P A R A M 3 , P A R A M 3 ) = H O T I I F ( M M A S S . G T . 1 ) J A C O B ( P ARAM 3 , P A R A M 3 ) = H D T I + { C A M P ( N ^ A S S ) + H D T * ( B + C ) ) / 3 R I N E R T { N M A S S ) C * * * C A L C U L A T I O N OF E L E C T R I C PART OF J A C O B I A N - LOWER R I G H T S Q U A R E * * * 6 0 3 T H E = P P O L E S * T H E N E W ( N M A S S ) W N E W D T = T H E - ( H O T * W C L C ( N M A S S ) * T H E O L D ( N M A S S ) ) * P P O L E S T H E 2 = 2 . D 0 * T H E A C 0 N 4 = A C C N 2 * 0 C 0 S ( T H E 2 ) , A C 0 N 5 = A C C N 2 * 0 S I N ( T H E 2 ) A C 0 M 6 = A C C N 3 * D C 0 S ( T H E 2 ) A C 0 N 7 = A C C N 3 * 0 S I N ( T H E 2 ) A C 0 N 8 = D C G S ( T H E ) A C C M 9 = D C C S ( T H F - A N G ) A C O N 1 0 = D C C S ( T H E + A N G ) A C 0 N 1 1 = D S I N ( T H E ) A C C N 1 ? = D S I N ( T H E - A N G ) A C 0 N 1 3 = D S I N ( T H E + A N G ) A C C N S 1 = A C C N 4 * A C G M A 4 - A C 0 N 7 * A C 0 N I B A C 0 N S 2 = A C C N 5 * A C C M A - A C 0 N 6 * A C C N l B B C 0 N S 1 = A C C N 2 * A C 0 M A - 0 . 5 D 0 * ( V Z B ( I 2 J - V Z R C I 3 ) - V Z C ( I 2 ) + V Z C ( I 3 ) ) B C O N S 2 = - A C O N 3 * ( V Z A ( I 2 ) - V Z A U 3 ) - V Z B ( U > * V Z 3 ( I 3 >+VZCf I 1 ) - V Z C ( 1 2 ) ) C 0 N I 0 A = A C C N 8 * C C N S A1 + A C 0 N 9 * C 0 N S A 2 + A C O N 1 C * C O N S A 3 CONI OB = ACCN ' l l * C G N S A l * A C O N 1 2 * C O N S A 2 * A C C N 1 3 * C C N S A3 D C O N S l = A C Q N * ( A C O N 8 * C 0 N S A 4 + A C C N 9 * C O N S A 5 + A C O N 1 0 * C O N S A 6 ) D C 0 N S 2 = A C C N * 1 A C O M 1* C 0 N S A 4 * A C 0 N 1 2 * CONS A 5 - A C 0 N 1 3 * C 0 N S A 6 » C 0 N S T A = C C N I 0 A * X ( 3 ) + A C O N l * ( A C 0 N 8 * E 1 ( 1 ) + A C O N 9 * E 1 { 2 ) + AC O N 1 0 * E 1 ( 3 ) ) V ( l ) = ( A C C N S 1 * B C 0 N S 1 ) * X 1 1 ) * ( A C O N S 2 * 8 C 0 N S 2 ) * X I 2 >+CONSTA C O N S T B = C O N T 0 B * X ( 3 ) • A C 0 N 1 * ( A C 0 N 1 1 * E 1 ( 1 ) + A C O N 1 2 * E 1 ( 2 ) + A C 0 N 1 3 * E l ( 3 ) V ( 2 ) = ( A C C N S 2 - B C O N S 2 ) * X ( 1 ) + ( - A C O N S 1 + R C O N S 1 1 * X ( 2 ) + C 0 N S T R V ( 3 ) = D C 0 N S 1 * X ( 1 ) - D C 0 N S 2 * X ( 2 ) * C 0 N I 0 C * X ( 3 ) * 1 . 0 0 / C O N S T 1* ( E l ( D + E K 2 3 4 E K 3 ) ) V ( 4 ) = V 0 L D ( 4 ) DO 6 C 5 J = 5 , M U 1 6 0 5 V ( J ) = O . D O N T E R M I N A L S Y S T E M F O R T R A N G ( 4 1 3 3 6 ) GEN 72. C * * * E V A L U A T E P R E S E N T F L U X L I N K A G E S * * * OO 6 0 7 J = l , * U l P S I ( J ) = O . D O DO 6 0 6 K=1,MU1 6 C 6 P S I ! J ) = P S I ( J ) - > - L P A R K { J , K > * X ( K I 6 0 7 C O N T I N U E TE = ( P S I ( 1 ) * X < 2 ) - P S I ( 2 ) * X ! 1 ) ) * P P O L E S B=O.DO C=O.DO I F ! I C O N V . E G . 1 ) GC TO 6 1 5 K = P A R A M 3 - 1 C * * * B E G I M F I L L I N G E N T R I E S O F E L E C T R I C P A R T I T I ME D E P E N D E N T ) * * * DO 6 C 8 I = N M A S S , K K 1 = I - P A R A M 1 DO 6 0 8 J = N M A S S , K K 2 = J - P A R A M 1 6 C 8 J A C O B ! I , J ) = L P A R K ! K 1 , K 2 J J1 = NMASS J2=NMASS+1 J3=NMASS + 2 K 1 = P A R A M 3 - 1 DO 6 1 0 J = N M A S S , K 1 K 2 = J - P A R A M l J A C O B ( J 1 , J ) = J A C O B ( J l , J ) * W N E W D T * L P A R K I 2 , K 2 ) J A C O R ( J 2» J ) - J A C O B ( J 2 , J J - W N F W O T * L P A R K ( 1 , K 2 ) 6 1 0 J A C O B ( J » J ) = J A C O R ( J , J ) + H O T * P V E C ( K 2 \ J A C O B ! J 1 , J 1 ) = J A C C B ( J l , J 1 ) + H D T * ( A C O N S 1 + RCONS 1 ) J A C O B ! J l , J 2 ) = J A C O B ( J l , J 2 J + H D T * ! A C 0 N S 2 + B C 0 N S 2 ) J A C O B t J l , J3 ) = H D T * C C N I O A J A C O B ! J l , P A R A M 3 ) = ! ( - 2 . 0 0 * A C 0 N S 2 * X { 1 ) + 2 . D 0 * A C C N S 1 * X ( 2 ) - C O N S T B ) * 3 HOT +P S I ! 2 ) ) * P P O L E S J A C O B I J 2 , J 1 ) = J A C C B ! J 2 , J l I * H D T * I A C 0 N S 2 - R C O N S 2 ) J A C O B ! J 2 , J 2 ) = J A C O B ( J 2 , J 2 ) + H D T * ! - A C C N S 1 + B C C N S I ) J A C O B ! J 2 , J 3 ) = H D T * C C N 10B J A C O B ! J 2 , P A R A M 3 ) = ( { 2 . 0 0* AC ON S 1 * X ! 1 ) + 2 . D O * A C C N S 2 * X { 2 ) + C O N S T A ) * 3 H O T - P S I !1 ) ) * P P O L E S J A C O B ! J 3 , J l > = H 0 T * C C C N S 1 J A C O B ! J 3 , J 2 ) = H 0 T * D C 0 N S 2 J A C O B ! J 3 , J 3 ) = J A C C B ( J 3 , J 3 ) + HD T *C ON I OC J A C O B ! J 3 , P A R A M 3 ) = ( - D C 0 N S 2 * X ! 1 ) + D C 0 N S 1 * X 1 2 ) ) * H D T * P P O L ES J = N M A S S + 3 DO 611 I = J , K 1 611 J A C O B ! ! , P A R A M 3 ) = C . D 0 C * * * C A L C U L A T E T I M E D E P E N D E N T ROW UNDER E L E C T R I C B L O C K OF J A C O B I A N * * * I F ( N M A S S . G T . l ) GC T O 6 0 9 R I N E R T l 1 ) = 1 . 0 1 0 D A M P ( 1 ) = 0 .DO 6 0 9 0 0 6 1 2 1=1,MU1 6 1 2 J A C 0 B ! P A R A M 3 , P A R A M 1 + I ) = L P A R K ! 1 , I ) * X ! 2 ) - L P A R K ! 2 , 1 ) * X ! I ) J A C O B I P A R A M 3 . J I ) = J A C O B I PARAM 3 , J D - P S I 1 2 ) J A C O B ! P A R A M 3 , J 2 ) = J A C O B ( P A R A M 3 , J 2 ) + PS I ! 1 ) DO 6 1 4 1=1,MU1 61A J A C O B I P A R A " 3 ,PAR AM I * I ) = J A C O B ( P A R A M 3 , P A R A M l ^ I > * H D T * P P O L E S / R I N E R T 3!NMA SS ) . C B E G I N C A L C U L A T I O N OF R I G H T HAND S I D E V E C T O R • G « * * * * * C * * * M E C H A N I C A L PART F I R S T * * * 6 1 5 I F I N R O T O R . E O . l ) GO T O 6 1 7 DO 6 1 6 I = 1 , P A R A M 2 N T E R M I N A L S Y S T E M F O R T R A N G ( 4 1 3 3 6 ) GEN 73. 6 1 6 T H E T A S I I ) = THENEW( I) 6 1 7 T H E T A S ( N R O T O R ) = T H E NEW(NMASS ) I F ( N ^ A S S . E O . 1 ) GO TO 6 2 2 I F ( N R O T O R . E C . N M A S S ) GO TO 6 1 9 K=NR .OTOP - l DO 6 1 8 I = K , N M A S S 6 1 8 T H E T A S ( I } = T H E N E W ( 1 - 1 ) 6 1 9 D = T H P T A S ( 1 ) DO 6 2 0 I=1,PARAM1 E = ST I F F ( I ) F = T H E T A S ( 1 + 1 ) HNEW2( I ) = - 8 * C + ( B - * E ) * 0 - E * F B = E C=D 6 2 0 D=F H N E W 2 1 N M A 5 S ) = - B * C + B * D B=HNEW2(NROTOR ) I F ( N R O T O R . E C . N M A S S ) GO TO 6 2 2 DO 6 2 1 J = N R 0 T 0 R , P A R A M 1 6 2 1 H N E W 2 ( J ) = H N E W ? ( J + l ) 6 2 2 H N E W 2 ( N M A S S ) = B * T E DO 6 2 3 J = l , N R O T O R 6 2 3 HNEW2 ( J ) = H N E W 2 ( J ) - T O R Q ' J E ( J ) DO 6 2 4 1 = 1 , N M A S S HNEW2( I ) = HNEW2I I ) * H D T / R I N E R T (I ) B = ( H D T I + C A M P { I l / R I N E R T t I ) ) * T r E N E W ( I ) HNEW1 ( I )=B C = B - H I S T 1 { I ) 6 2 4 G ( I ) = - C - H N F W ? ( I ) - H I S T 2 ( I ) G ( P A R A M 3 ) = G ( N M A S S ) I F ( I C O N V . E Q . l ) GO TO 6 4 2 C B E G I N R I G H T HAND S I D E OF E L E C T R I C A L E Q U A T I O N S * * * * * L = P A R A M 3 - 1 DO 6 2 5 J = N M A S S , L K = J - P A R A M 1 6 2 5 G ( J ) - - P S n K ) + P S I O L D ( K ) - H D T * ( V ( K ) + V O L D ( K ) + R V E C ( K ) * ( X ( K ) + X O L D ( K ) ) ) C = T H E N E W ( J 1 ) - T H E C l C ( J l > G t J l ) = G ( J l ) - P S I ( 2 ) * P P 0 L E S * C + H D T * P P 0 L E S * W 0 L D ( J l ) * ( P S I ( 2 ) - P S I O L D ( 2 : G ( J 2 )=G( J 2 ) + PS K 1 ) * P P O L E S * C - H D T * P P O L E S * W O L D ( J 1 ) * ( P SI ( 1 ) - P S IOL D ( 11 C * * * S O L V E T H E S E T OF L I N E A R E Q U A T I C N S ( U P T O 15 X 1 5 ) * * * C E T = 1 . D - 1 4 I J K = 0 C A L L D S L I M P ( J A C 0 B , C A , G , D E L T A X , R 7 , I P E R M , P A R A M 3 . 1 5 , D E T , 1 , 1 4 ) I F ( D E T . L T . O . D O ) GO TO 631 DO 6 3 3 1 = 1 , P A R A M 3 6 3 3 I F ( D A B S ( D E L T A X ( I ) ) . G T . l . D - 9 ) I J K = I J K + i I F ( N M A S S . E Q . l ) GC TO 6 3 7 DO 6 3 8 J = 1 , P A R A M 1 6 3 8 THENEVM J ) = T H E N E W ( J ) + D E L T A X ( J ) 6 3 7 T H E N E W ( N M A S S ) = T H E N E W ( N M A S S ) +DEL T A X ( P A R A M 3 ) , D O 6 3 9 J = 1 , M U 1 K 1 = P A R A M 1 + J 6 3 9 X( J ) =X( J ) + D E L T A X ( K l ) I TER = I T E R + 1 I F ( I T E R . G T . 1 0 ) W R I T E ( 6 , 6 5 6 ) 6 5 6 F O R M A T { • • , » * * * * * S O L U T I C N NOT W I T H I N E P S I L O N A F T E R 10 I T E R A T I O N S . 3 L A S T V A L U E S U S E D * * * * * ' ) I F ( I T E R . G T . 1 0 ) GO TO 6 4 1 N T E R M I N A L S Y S T E M F O R T R A N 0 ( 4 1 3 3 6 ) G E N I F ( I J K . G T . O ) GO TO 6 0 3 W R ! T E ( 6 , 6 4 0 I I T E R 6 4 0 FORM AT ( • ' , • S O L U T I O N BY NEWTCN* ' S METHOD C O N V E R G E D IN ' , I 2 , ( IT E f (DTICNS • ) 641 I C 0 N V = 1 GO T O 6 0 3 6 4 2 DO 6 4 5 1=1 ,NMASS WOLD ( I ) = H D T I * T H E N E W ( I ) - ( W O L O ( I ) + H D T I * T H E O L C ( I ) ) T H E O L O ( I ) = THENEW( I ) HI ST 1 ( I ) = H N E W 1 ( I t - 2 • D O * W C L D ( I ) 6 4 5 H ! S T 2 ( I ) = H N E W 2 ( I ) 1 L M = ( N - l ) *7 DO 6 4 4 J = 1 , M U 1 Z T R ( I L M + J ) = X ( J ) Z T X ( I LM + J ) = V ( J ) 6 4 4 Z T Y ( I L M + J ) = P S I ( J ) C NOW C O N V E R T V O L T A G E S AND C U R R E N T S TO P H A S E DOMAIN W = A C 0 \ ' 1 * A C C N 1 4 * X ( 3 ) A C 0 N 1 A = A C C N 1 * ( A C C K 8 * X ( I) + A C O M 1*X ( 2 ) ) + W A C 0 N 1 B = A C C N 1 * ( A C C N 9 * X ( 1 ) + A C 0 N 1 2 * X ( 2 ) ) + W A C 0 N S 1 = A C C N 1 * ( A C C M O * X ( l >~AC0N13*X (2 ) )*W W = A C 0 N 1 * « C 0 N 1 4 * V ( 3 ) A C 0 N 2 A = A C C N 1 * ( A C C N 8 * V ( 1 ) + A C 0 N 1 1 * V ( 2 >)*W A C 0 N 2 B = A C C N l * ( A C 0 N 9 * V ( 1) + A C 0 N 1 2 * V( 2 ) ) + W A C 0 N S 2 = A C C N 1 * < A C C N 1 0 * V ( 1 ) + A C 0 N 1 3 * V { 2 ) ) +W TE = P P Q L E S * ( P S I ( 1 > * X ( 2 ) - P S I (2 ) * X ( 1 ) I C U R R T ( 1 ) = - A C O N l A C U R R T ( 2 ) = - A C O N l B C U R R T ( 3 ) = - A C C N S l W R I T E ( 6 , 6 5 0 ) A C 0 N 1 A , X ( 1 ) , X ( 3 ) , X ( 4 ) , T E , W O L D ( N M A S S ) 65 0 FORMAT ( • • , M - A = • , CI 4 . 7 , 2X , « 1-0= • , DI 4 . 7 , 2X , • I - 0= » , DI 4 . 7 , 2X , « I -F = a D 1 4 . 7 , 2 X , « T - E = , , 0 1 4 . 7 , 2 X , « R 0 T 0 R S P E E D = « , 0 1 4 . 7 ) i T H E T ( N) = P P O L E S * T r - ENEW ( NM ASS ) - O M E G A * T I M E - 1 . 5 7 0 7 9 6 3 3 D 0 B=O.DO C=O.DO 0=0.DO I F ( N M A S S . F O . l ) GO T O 6 5 1 R= ( T H E T A S ( P A R A M 2 ) - T H E T A S ( N R O T O R ) ) * S T I F F { P A R A M 2 ) I F ( N M A S S . E 0 . 2 ) GC T O 6 5 1 C = ( T H E T A S ( N P 0 T 0 R - 2 ) - T H E T A S ( P A R A V 2 ) )* ST I F F < N R Q T O R - 2 ) I F ( N M A S S . G T . N R O T O R ) D = ( T H E T A S ( N R O T O R I - T H E T A S ( N R O T O R + 1 ) ) * S T I F F ( N R C aoR) 6 5 1 W R I T E < 8 ) A C 0 N 1 A , A C O N 1 B , A C O N S 1 • A C 0 N 2 A , A C 0 N 2 P , A C C N S 2 , a X ( 4 ) , T E , B , C , D , T H E T ( N ) GO TO 6 5 5 6 3 1 W R I T E ( 6 , 6 5 3 ) 6 5 3 F O R M A T ! • « , ' S O L U T I C N F A I L E D , C H E C K YOUR M A T R I C E S ' ) 6 5 5 C O N T I N U E I = I S A V E RETURN END riCNS IN E F F E C T * 1 0 , E 8 C D I C , S O U R C E , N O L 1 S T , N O D E C K , L O A D , N O M A P ("IONS IN E F F E C T * NAME = GEN , L I NEC NT = 6 0 4 T I S T I C S * S O U R C E S T A T E M E N T S = 2 4 9 , PROGRAM S I Z E = 1 3 0 3 2 V T T S T I C S * NO D I A G N O S T I C S G E N E R A T E D *S IN GEN 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065483/manifest

Comment

Related Items