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.

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 p r e s e n t i n g t h i s  thesis  in p a r t i a l  f u l f i l m e n t o f the requirements f o r  an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, the L i b r a r y s h a l l make i t I  freely available  f u r t h e r agree t h a t p e r m i s s i o n  for  I agree  r e f e r e n c e and  f o r e x t e n s i v e copying o f  this  that  study. thesis  f o r s c h o l a r l y purposes may be granted by the Head of my Department or by h i s of  this  representatives. thesis  It  i s understood that copying or p u b l i c a t i o n  f o r f i n a n c i a l gain shall  written permission.  Department of  pjftrtirica  The U n i v e r s i t y o f B r i t i s h  2075 Wesbrook Place Vancouver, Canada V6T.1W5  Date  May  8,  13 7 8  Columbia  not be allowed without my  ii  ABSTRACT  In  recent years,  interaction  in  interaction  between  in  by  turbine-generator  transmission  these  lines,  disturbances  turbine-generator to  represent  than  as  shaft.  voltage  Addition  trical  is  T h e model in  shaft  the machine.  phase  plied the  to  cases  conclusions  a study were  in  system  are  these  system series  Shaft  compen-  destroy  t h e r e was  programs  occurs  torques  damage o r  cases,  model  to  and  a  possesses  large  and no m a g n e t i c Program  the e x t e r n a l rule  followed  use  by  both  saturation is  through  electric  implicit  the  d e g r e e of  detail  a  the  to  verify  as  to  method  t h e new s o l u t i o n  seen is  the synchronous  synchronizations  of  an u n l o a d e d  the  findings  of  machine  generator,  a manufacturer.  ap-  solve  technique,  of  iron  three-  integration  of N e w t o n ' s  elec-  generality.  of  network  applications  agreement w i t h  the  Transients  equations.  presented  caused  a need  i n more  the U.B.C.  represents  Trapezoidal  of p o s s i b l e  to  transients  the T r a n s i e n t s  algebraic  as  in  the g e n e r a t o r  c i r c u i t of  of f a u l t y  severe  study  torsional  reactances.  excitation to  may a r i s e  synchronization.  generator  of  This  and e l e c t r i c a l  that  computer  dynamics  terminals.  nonlinear  to  behind  T h e model  Interface  by e x a m p l e s In  faulty  digital  the machine e q u a t i o n s  Test followed  in  sources  Thevenin equivalent  resulting  model.  in  constant  from the g e n e r a t o r  power s y s t e m s .  phenomenon,  Therefore,  torsional  assumes  in  shaft  of a s y n c h r o n o u s  described.  and  and  arisen with electro-mechanical  may be s u f f i c i e n t l y  generators  simple  Program  resonance  have  units  the mechanical  the subsynchronous  sated  problems  i ii  TABLE OF CONTENTS  Page ABSTRACT  ii  LIST OF ILLUSTRATIONS  v  ACKNOWLEDGEMENT  vi  1 . INTRODUCTION. 2.  SYNCHRONOUS  .  1  MACHINE MODELLING  . .  k  2.1  M o d e l l i n g of E l e c t r i c a l E f f e c t s  k  2.2  I n t e r f a c e to the E l e c t r i c a l Network  11  2 . 3 M o d e l l i n g o f Mechanical E f f e c t s 3.  STEADY STATE 3.1  4.  12  INITIALIZATION OF MACHINE EQUATIONS.  . . .  E l e c t r i c a l Variables  16 16  3 . 2 Mechanical System V a r i a b l e s  18  THE TIME STEP - NUMERICAL SOLUTION OF THE MACHINE EQUATIONS  20  4.1  Choice of  Integration  Method  20  4 . 2 S o l u t i o n o f the Machine Equations by Newton's Method. 2 0 4.3 5.  23  STEADY STATE AND TIME-STEP SOLUTION 5.1  6.  Form of the J a c o b i a n PROGRAMS  P r e l i m i n a r y Remarks About the Programs  .31 31  5 . 2 Steady S t a t e - S u b r o u t i n e "MSTEA"  33  5 . 3 Time Step - S u b r o u t i n e  3^  ,:I  GEN"  5 . 4 Flow Diagram o f S o l u t i o n by Newton's Method  36  TESTING OF THE SIMULATION METHOD  38  6.1  Initial  6 . 2 Single 6.3  Tests  38  Line-to-Ground Fault  IEEE Subsynchronous  Resonance Benchmark T e s t Case  39 . . 40  iv  7.  EXAMPLES OF PROGRAM APPLICATIONS 7.1  Sustained Generator  7.2 8.  J  43  Three-Phase F a u l t f o r an Unloaded ..  . . •  F a u l t y S y n c h r o n i z a t i o n o f an Unloaded Generator  43 . . .  CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH.  . .  43 48  REFERENCES  50  APPENDIX I  52  APPENDIX II  63  APPENDIX III  69  LIST OF ILLUSTRATIONS Fi gure  Page  1  E l e c t r i c a l model of synchronous  2  Terminal c o n n e c t i o n of the network to the machine.  12  3  Mechanical  14  4  Model 1ing o f mechanical i-th  machine.  7  r e p r e s e n t a t i o n of a t u r b o g e n e r a t o r .  rotational  system  in the v i c i n i t y o f  the  mass.  14  5  Phasor diagram o f the synchronous  machine  in steady  state.  6  Grouping o f equations  7  Nonzero s t r u c t u r e o f the J a c o b i a n m a t r i x f o r an eight-mass seven-winding generator model (producing the l a r g e s t matrix)  30  8  Flow diagram o f s o l u t i o n  37  9  Simple network c o n f i g u r a t i o n .  and independent v a r i a b l e s .  by Newton's  25  method.  38  .10  Network c o n f i g u r a t i o n  for single  11  E l e c t r i c network f o r the s i m u l a t i o n resonance. model o f the s h a f t  17  line-to-ground of  fault.  39  subsynchronous 40  12  Mechanical  system.  13  Subsynchronous resonance, phase c c u r r e n t ( u p p e r ) , E l e c t r o magnetic torque 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 torque ( l o w e r ) .  42  14  System diagram f o r t h r e e - p h a s e f a u l t .  44  15  E l e c t r o m a g n e t i c torque on r o t o r mechanical torque ( l o w e r ) .  41  ( u p p e r ) , and  LPB-rotor 44  16  System diagram f o r f a u l t y s y n c h r o n i z a t i o n .  46  17  Torques due to f a u l t y s y n c h r o n i z a t i o n  (3= 120°)  47  18  Torques due to f a u l t y  (3= 240°)  47  synchronization  vi  ACKNOWLEDGEMENT  I would l i k e to extend my s i n c e r e thanks  to my s u p e r v i s o r s ,  Dr. P. Fenwick and Dr. H. Dommel f o r t h e i r guidance, a d v i c e and u n r e m i t t i n g encouragement and p a t i e n c e d u r i n g my p e r i o d of graduate s t u d y .  I a l s o wish  to express my deep a p p r e c i a t i o n to Dr. V. Brandwajn f o r h i s enthusiasm this  p r o j e c t , his  used in t h i s  t i m e l y a d v i c e , and f o r s u p p l y i n g some of the t e s t  for  data  thesis.  The f i n a n c i a l  support of the N a t i o n a l  and the U n i v e r s i t y o f B r i t i s h Columbia  Research C o u n c i l o f  is g r a t e f u l l y  acknowledged.  Canada  1.  1.  This  thesis  serves  INTRODUCTION  t h r e e main purposes.  The f i r s t  develop and then to add the e l e c t r o - m e c h a n i c a l synchronous  purpose  in the form of an a c c u r a t e and reasonably e f f i c i e n t s i m u l a t i o n is  to show t h a t  those of o t h e r methods. that  results  The f i n a l  computer programs  Program  routine.  o f the new program are  in agreement w i t h  purpose o f  is  this  thesis  to demonstrate  the program may be used to study c e r t a i n problems of c u r r e n t Digital  to  machine dynamics  to the framework o f the e x i s t i n g U.B.C. 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  The second  is  s p e c i f i c a l l y designed  f o r the  importance. simulation  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 phenomena in e l e c t r i c networks, have become a very powerful of  general-purpose  tool  in power system a n a l y s i s .  the power system can be r e p r e s e n t e d w i t h a high degree of  Some elements sophistication,  Others are o v e r s i m p l i f i e d , e i t h e r because more d e t a i l e d r e p r e s e n t a t i o n  is  not warranted f o r the p a r t i c u l a r case s t u d y , o r , because c e r t a i n n o n l i n e a r effects  are too c o m p l i c a t e d to model.  practice  to model a synchronous  grams as  sinusoidal  R  s  + jXV , w i t h R d ' s J  transient  voltage  machine  sources  r e c e n t l y i t had been the  in 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  E|^  s a x  ' (<JJt+<j>) behind  rotor fluxes  subtransient  Some types o f p a r t i c u l a r l y true is  results  transients  persist  system.  subas-  The s i m p l e  d u r i n g a few c y c l e s a f t e r involving  fast  in-  t r a n s i e n t s [1].  f o r more than a few c y c l e s .  This  in cases where, f o l l o w i n g a d i s t u r b a n c e such as a f a u l t ,  i n t e r a c t i o n between the mechanical  electrical  s a l i e n c y and  a f t e r the t r a n s i e n t d i s t u r b a n c e .  i t i a t i o n of a disturbance for c e r t a i n studies  there  impedances  n  T h i s approach n e g l e c t s  machine model g i v e s q u i t e s a t i s f a c t o r y  is  pro-  b e i n g the armature r e s i s t a n c e and XV the d i r e c t - a x i s • d  reactance.  sumes c o n s t a n t  Until  Examples  t u r b i n e - g e n e r a t o r system and the  i n c l u d e the torque a m p l i f i c a t i o n caused by s u b -  synchronous nization. 1970  and  Nevada  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  In t h e wake o f t h e two s u b s y n c h r o n o u s 1971,  resulting  [2], u t i l i t i e s  i n s h a f t d e s t r u c t i o n a t t h e Mohave p l a n t  recognized  t h e need  spans  l o n g e r t h a n two o r t h r e e c y c l e s .  of  a program  Pacific  Gas  and  In 1974,  c o n t a i n i n g both e l e c t r i c a l  synchronous  machines.  T h e i r machine  resulted  and s h a f t  model was  v e r s i o n o f t h e U.B.C. p r o g r a m ) b u t has version[3].  T h i s machine  model was  recently  in  torsional  implemented  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  analytic  generators over time  a joint effort  Electricity  of  in southern  t o d e v e l o p more d e t a i l e d  to study the t r a n s i e n t b e h a v i o r of synchronous  E d i s o n and  ^synchro-  resonance, i n c i d e n t s  tools  California  faulty  program  by  Southern  development  dynamics with  of  Bonneville  (a much  expanded  been s u p e r s e d e d by a  faster  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  this  period [4], The  U.B.C. T r a n s i e n t s P r o g r a m was  commodate a s y n c h r o n o u s  g e n e r a t o r model,  b r a n c h w i t h n o n l i n e a r and terminal  bus  itialization  to ground. and  t e s t i n g o f the  The  major  thesis  and s h a f t  to ac-  that connects the generator  r e s e a r c h e f f o r t was  the t i m e - s t e p s o l u t i o n ,  proceeds w i t h development  torsional  s o l u t i o n s w h i c h was  dynamics  the machine  F o r t r a n code  directed  t o the i n -  i m p l e m e n t a t i o n , and  t o compare w i t h  machine  under c o n t r o l  new  b a s e d on N e w t o n ' s m e t h o d .  i s shown f o r t h e s t e a d y s t a t e a n d  implemented  then d e s c r i b e s examples  o f the synchronous  f o l l o w e d by a p p l i c a t i o n o f a  e q u a t i o n s [6]  Test cases are presented to e s t a b l i s h t e c h n i q u e and  [5]  as a t h r e e - p h a s e n e t w o r k  time-varying parameters  particularly  technique f o r s o l v i n g The a c t u a l  treated  modified  program.  This electrical  previously  the time-step  of the Transients  the r e l i a b i l i t y  o f t h e new  r e s u l t s o b t a i n e d by o t h e r m e t h o d s .  where a p p l i c a t i o n o f the program  i s used  Program. solution The  thesis  in studying  3.  problems such as  the s h a f t  torques caused by a three-phase f a u l t and f a u l t y  s y n c h r o n i z a t i o n of an unloaded g e n e r a t o r . mendations are given f o r f u r t h e r r e s e a r c h .  F i n a l l y , conclusions  and recom-  4.  2. SYNCHRONOUS MACHINE MODELLING  2.] M o d e l l i n g of E l e c t r i c a l A mathematical obtain  Effects  model of a p h y s i c a l  i n f o r m a t i o n about performance and c h a r a c t e r i s t i c s without having  t e s t the a c t u a l d e v i c e in the l a b o r a t o r y o r ematical  o f the d e v i c e .  The u s e f u l n e s s  is apparent when one c o n s i d e r s s e r v i c e a large  and p h y s i c a l  impossibility  available.  behavior  device modelling in removing  its  device involves  The model s t r u c t u r e  parameters of the model  is of no use to s e l e c t a h i g h l y  a certain  is d e r i v e d from  always the  possibility  to whatever data  sinusoidal practical  voltage  For example, one could not study  to employ a machine model  e f f e c t s as eddy c u r r e n t s  the change  s h o r t c i r c u i t w i t h a generator  sources behind r e a c t a n c e s .  Similarly,  consistent  representation of i t would be imtransient  in the r o t o r when the time p e r i o d of study in the o v e r a l l  from  in f i e l d c u r -  r e p r e s e n t i n g such c o m p l i c a t e d  that such phenomena play no r o l e  is  d e t a i l e d and c o m p l i c a t e d  However, the degree o f c o m p l e x i t y of the model must be  purpose.  from  purposes.  understanding but t h e r e is  of f i t t i n g  rent f o l l o w i n g a t e r m i n a l  short  reality.  ( l i n e a r o r non-  i f the parameters cannot be c o m p l e t e l y determined o r e s t i m a t e d  test data. with  It  the d i f f i c u l t y and expense  to  The term " m a t h -  the b a s i c p h y s i c a l  and convenience of p h y s i c a l  model of a p h y s i c a l  degree of s i m p l i f i c a t i o n of observations  that d e s c r i b e s  turbogenerator for test  Any mathematical  model  in the f i e l d .  model" i m p l i e s an e q u a t i o n or system o f e q u a t i o n s  l i n e a r , d i f f e r e n t i a l or algebraic)  or  d e v i c e is o f t e n employed to  is  so  behavior of the p h y s i c a l  system. To be o f p r a c t i c a l  a p p l i c a t i o n , a model of a synchronous  should accommodate two kinds of  rotor:  generator  5.  (1) S a l i e n t  pole r o t o r used in slow speed h y d r o - e l e c t r i c g e n e r a t i n g  plants; (2) s o l i d l y - f o r g e d  steel  round r o t o r  in uniform a i r - g a p  turbogenerator  un i t s . Salient  pole synchronous  g e n e r a t o r s o f t e n have more than one pole p a i r ;  lower the head o f water behind the dam, the l a r g e r Also,  the number o f  these machines f r e q u e n t l y have damper or " a m o r t i s s e u r "  These damper c i r c u i t s , t a k i n g  switching.  poles.  circuits.  the form of brass o r copper bars housed in the  laminated p o l e f a c e s , serve s e v e r a l ing of r o t o r o s c i l l a t i o n s  the  useful  purposes  following disturbances  [7], i n c l u d i n g the damp-  such as s h o r t c i r c u i t s or  Uniform a i r - g a p machines, on the o t h e r hand, are u s u a l l y  and s l o t s are m i l l e d i n t o the s o l i d  two pole  round r o t o r body to support the f i e l d  winding. With the above o b s e r v a t i o n s , chronous machine w i l l assumptions  the mathematical model of the s y n -  be s e t f o r t h w i t h i n the framework o f the  following  [5]:  (1) A l l  f l u x d e n s i t i e s are p r o p o r t i o n a l  so that  to the c u r r e n t s producing them  there is no s a t u r a t i o n o r o t h e r n o n l i n e a r e f f e c t s .  (2) The machine has a b a l a n c e d , smooth s t a t o r and the s t a t o r windings are Y-connected with the n e u t r a l p o i n t grounded. (3)  Eddy c u r r e n t paths  in the r o t o r body o f u n i f o r m a i r - g a p machines o r  in the a m o r t i s s e u r c i r c u i t s of s a l i e n t  pole machines can be r e p r e -  sented by c o n c e n t r a t e d h y p o t h e t i c a l w i n d i n g s , axis  (D) , the o t h e r  thetical  winding  ing eddy c u r r e n t s  (g)  in the quadrature a x i s in the quadrature a x i s  one in the d i r e c t  (Q_) , and another hypoto r e p r e s e n t deep f l o w -  in the case o f round r o t o r machines.  (4) Each s e l f and mutual  inductance may be expressed as a sum of a  6.  c o n s t a n t and a s i n u s o i d a l 0 is  20(t). the  f u n c t i o n o f r o t o r angle 0(t)  the angle between phase " a " and the d i r e c t a x i s  f i e l d winding measured  (5) Mutual  and/or  in the d i r e c t i o n o f  of  rotation.  inductances between armature and f i e l d and between armature  and damper windings are e q u a l . (6) During  the time of the s i m u l a t i o n ,  steady s t a t e  f i e l d voltage  i n t o the e x t e r n a l  c u r r e n t flow  network.  carrying different currents.  pole machines c o n s i s t s  Exact  r e p r e s e n t a t i o n would  q u i r e a l a r g e number o f c o n c e n t r a t e d h y p o t h e t i c a l is  sufficient  coils  on the quadrature a x i s  the eddy c u r r e n t paths  winding having an i n f i n i t e number o f  1  f  1  f  1  , g , D ,and 1  1  1  1  'Q_'.  1  a , b , c 1  All  1  in phase  1  1  1  of t h r e e i d e n t i c a l ,  and f o u r unequal  machine  is  sources.  shown in F i g .  A 1.  coordinates.  -[R]T - ^ ?  = -[R]T the v e c t o r o f f l u x  nous machine  lumped  laws may be a p p l i e d , y i e l d i n g one e q u a t i o n per  7=  ijj is  synchro-  r o t o r windings except the f i e l d  schematic r e p r e s e n t a t i o n o f the synchronous  coil  purposes  bars.  are s h o r t c i r c u i t e d as they are not connected to v o l t a g e  K i r c h ' h o f f ' s and F a r a d a y ' s  re-  in the iron can be c o n s i d e r e d as a cage  p l a c e d armature windings  r o t o r windings  of  on the d i r e c t  For round r o t o r  [8],  The model o f the e l e c t r i c a l system c o n s i s t s symmetrically  but f o r most  to represent the damper winding as one c o i l  a x i s and another c o i l nous machines  its  is out from the g e n e r a t o r and  The c a g e - l i k e damper winding of s a l i e n t  it  constant at  value.  (7) The d i r e c t i o n of s t a t o r  many c i r c u i t s  is  linkages  in phase c o o r d i n a t e s  [8].  (l)  {[L]T} The mathematical  (2) model o f the s y n c h r o -  is d e f i n e d by the seven f i r s t o r d e r  dif-  7.  Fig.  1.  E l e c t r i c a l model of synchronous  machine.  8.  ferential  equations  represented by (l)  by o t h e r s and r e s u l t s  or  (2).  T h i s model has been t e s t e d b  agree q u i t e w e l l w i t h f i e l d t e s t s  [k],  be adequate f o r most types 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 In d e r i v i n g  (1) or  [3]  and should  s t u d i e s [10].  (2), the s t a t o r equations were expressed  frame of r e f e r e n c e a t t a c h e d to the s t a t o r and the r o t o r equations a r a t e frame moving with the r o t o r .  that are f u n c t i o n s of r o t o r a n g l e .  i n d u c t a n c e s , c o n t a i n components  Owing to the r o t o r - a n g l e dependence of  the inductance m a t r i x , the d i f f e r e n t i a l equations dependent parameters, making s o l u t i o n by using o n l y one frame of  thereby a v o i d i n g of v a r i a b l e s formation  the time^varying  that accomplishes  [11].  It  in a s e p -  The s e l f - i n d u c t a n c e s o f the s t a t o r wind-  ings t o g e t h e r w i t h the i n t e r - w i n d i n g mutual  results  in a  given by (2) c o n t a i n t i m e -  r a t h e r awkward.  A great  r e f e r e n c e and a t t a c h i n g inductance c o e f f i c i e n t s .  this  change  is u s u a l l y  d e f i n e s a new s e t o f s t a t o r  this The  simplification to the r o t o r , transformation  called Park's  voltages,  trans-  c u r r e n t s , and  f l u x e s from a p r o j e c t i o n o f the a c t u a l winding v a r i a b l e s on t h r e e a x e s ; one along  the d i r e c t a x i s of the f i e l d w i n d i n g ,  c a l l e d the d i r e c t or " d "  axis,  a second along an a x i s p e r p e n d i c u l a r to the f i e l d w i n d i n g ,  c a l l e d the quadra-  ture or " q " a x i s ,  The t r a n s f o r m a t i o n  and the t h i r d on a s t a t i o n a r y  of s t a t o r c u r r e n t s and v o l t a g e s v  "0" a x i s .  to the new c o o r d i n a t e system is d e f i n e d by: p  P where [T]  (3)  = [T]v [T]i P  1  0  (5)  o; i COS0  cos (0-120°)  cos(0+]20°)  (6) /2  /2  /2  9.  which  is  slightly  as b e f o r e , [P]  is  is  d i f f e r e n t to P a r k ' s o r i g i n a l  transformation.  the angle between phase " a " and the d i r e c t a x i s .  Angle  6,  The m a t r i x  orthogonal: [Pf  and the t r a n s f o r m a t i o n  is  P  (7)  1  Rewriting  [12].  [R]T  the t r a n s f o r m a t i o n v  [P]  power i n v a r i a n t  v = and a p p l y i n g  =  1  = - [L  (2) as  (8) gives  At -  P dt p J  [5]: [R  ]Tp  -  [L ] i p p  (9)  1  where •'! L  c  0  L  0  0  0  0  L  c  0  0  c  0  J  0  fa  fa  fa  0  0  L  f  M  f  fa  0  0  M  f  L  D  °  ^  °  J"q  0  0  0  0  0  0 {[Pf  dt'  0  L  0  0  q  0 -fa  d = 10  0  0  0  0  0  0  o  fa  0  42  q 42 q  0  0  0  0  0  0  (JO)  M  q  L  M q L  P  V2 q  9  ]> \I2 q  0  0  0  0  (11)  10.  or = -R  d  v  q  = -R  o  v  =  dt*0 d .  dVf  " f  f'"  " D  D  , dt^D  - Q  o_-~  d , dt^Q  R  q  dt^q d .  o"  s  f  (12)  d  0  =  0  =  0  = -R  R  R  dt g v  g g  i s d e f i ned by:  i nkages, =  dt*d - u^  s d  = -R * s q  V  v  d .  i  L  d  !  d ^ +  ibt " f J +  +  L i  q q -.-\  *o -  W0  *f -/^f'd f'f f D /fVd V f D D ^0. fh\ VQ Vq + L  ij) r  g  Equations  useful  +  =  +  =/|M i  V2  (12) t o g e t h e r w i t h  nous machine model is  =  in P a r k ' s  q q  !  + H +  L  (13)  !  +  + M i + L i q Q gg  (13) d e s c r i b e the e l e c t r i c a l  part of  coordinates.  (12) and (13), i t  In a d d i t i o n to  to develop the e l e c t r o m a g n e t i c torque on the The  power at  the g e n e r a t o r t e r m i n a l s P  Substituting  p e  d  , + .i .  rotor.  is [13]:  = v ,i , + v I d d q q  from (12) f o r v . and v 'd - ~ q d  (1*0  gives:  Jt* ' *d <,--+ d>- s<  --"di*d  q  the synchro-  +w(  i  i  q  R  q  i d  + i q  )  (15)  .2 . 2 .2, Of the three groups o f  terms  resistance;  in  (15),  R s  in  the s t a t o r  of  change o f magnetic e n e r g y ; and t*>(^ i  ('  + d  (i , + i d dt d q at q d  q  ' ) q  is  - i|» i j )  r e p r e s e n t s the power a measure o f  the time  r e p r e s e n t s the power  loss rate trans-  11  ferred across the a i r gap from rotor to stator.  3  The a i r gap torque, T  »  r o t o r  can be defined as: T  The mechanical  e rotor  =  air-gap power mechanical rotor speed  and rotor speed, as™, are related to the e l e c t r i c  angle, 8 , m  angle, 0, and frequency, w, of the network by: 0 = n 9 P  (17-a)  m  co = n co P  (l7 b)  m  -  where n^ is the number of pole pairs in the generator.  Therefore, the elec-  tromagnetic torque i s : T 2.2  . = n (if/.i - ill i ,) rotor p d q - q d  (18)  6  T  T  Interface to the E l e c t r i c a l Network The seven f i r s t - o r d e r d i f f e r e n t i a l 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 e l e c t r i c a l 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 currents, allowing machine terminal voltages to be expressed 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 a V  b' c ' V  a S  in terms of stator  s  n  o  w  n  b  while the bus voltages are the stator voltages v , c  a  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 e f f e c t of the network on the machine is identical to the effect of  12.  a t h r e e - p h a s e Thevenin e q u i v a l e n t c i r c u i t connected to the machine  terminals,  Or: =  \M where the s u b s c r i p t  C "N^ N ^ N  +  J  e  No  (19)  ( t )  " N " denotes network q u a n t i t i e s  and " o " denotes open  c i r c u i t quanti t i e s .  Thev "N 7  0 0 e j 1111  i  11111  GENERATOR  Fig.  All  that  2. Terminal  remains  pressions  for v., d'  K  (12).  is  connection of  to t r a n s f o r m  i ,, i d'  0  coordinates, obtaining  ex-  , i_. which can be s u b s t i t u t e d q'  into  0  T h i s procedure is shown irn Appendix 1.  2.3 M o d e l l i n g  o f Mechanical  The synchronous rotating  Effects  machine mechanical  system c o n s i s t s  components coupled together by s h a f t s ,  the e l e c t r i c a l is  the network to the machine.  (19) to d-q-0  v , v„ in terms o f q'  c  ELECTRIC NETWORK  p a r t o f the g e n e r a t o r ,  as shown in F i g .  3-  of  As w i t h  the complexity o f the mechanical model  c o n t r o l l e d by the amount of a v a i l a b l e  the model.  of a s e r i e s  data and the intended purpose of  Many o f the damping c o n s t a n t s ,  f o r example, are very d i f f i c u l t  o b t a i n and are o f t e n not even'measured by m a n u f a c t u r e r s .  It would be q u i t e  to  13.  impractical  to develop a machine mechanical model w i t h ,  p h i s t i c a t e d nonviscous the parameters as  damping and n o n l i n e a r t o r s i o n a l  could be measured, such a d d i t i o n a l  for  instance,  springs.  complexity  Even i f  is  masses was s e l e c t e d as chronous machine.  the most general  The a c t u a l  number of  u l a t i o n depends on the i n d i v i d u a l pressure  (HP)  system having  mechanical  (IP)  loads  sim-  a high  t u r b i n e , three low (EXC), and a f e e d -  r o t a t i n g machinery.  However,  application  o f t h i s model goes beyond the r e p r e s e n t a t i o n of a t u r b o g e n e r a t o r . equally well  lumped  representation of a syn-  (LPA, LPB, LPC) , the r o t o r , an e x c i t e r  water pump (PUMP) o r o t h e r p i e c e o f  up to e i g h t  T h i s model accommodates  t u r b i n e , an i n t e r m e d i a t e p r e s s u r e  p r e s s u r e stages  behavior.  lumped masses used in a given  case.  all  unwarranted  i t does not c o n t r i b u t e to the understanding o f the b a s i c system A l i n e a r multimass-spring-damper  so-  It  could  r e p r e s e n t , f o r example, a motor d r i v i n g one or more mechanical  such as pumps, o r a combination o f a machine, t u r b i n e s , and mechanical  loads on a s h a f t . The f o l l o w i n g model of (1)  the mechanical  assumptions  are made in d e r i v i n g  system.  The system to be modelled c o n t a i n s  (2) Mutual is  damping  the mathematical  is  up to e i g h t  rotational  c o n s i d e r e d unimportant and only speed  masses. self-damping  used.  (3) The r e l a t i v e angular  deflections  in the s h a f t  sections  are  enough so t h a t the lumped masses may be connected by l i n e a r  small torsional  springs. (4) The p e r i o d o f s i m u l a t i o n  This  l i m i t a t i o n was  presentation system  is  is  short  imposed only because  s u f f i c i e n t l y general  l i k e l y to be e n c o u n t e r e d .  compared to the time constant  i t was  felt  that an eight-mass  to accommodate any a c t u a l  mechanical  of  re-  14.  HP  IP  LPA  LRB  LPC  Fig.  3- Mechanical  Fig.  4. M o d e l l i n g of mechanical system of the i - t h r o t a t i o n a l mass.  ROTOR  r e p r e s e n t a t i o n of a turbogenerator.  in the v i c i n i t y  . EXC  PUMP  15.  the governor system so that F i g . k shows the i - t h in r o t a t i o n a l  rotational  input mechanical  mass and i t s  2  i  2  d  vicinity.  form a p p l i e d to each lumped mass y i e l d s  order 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 .  d e. J.——-  torques are Newton's  second  a series of  second  i  M-1  D T  law  For the i - t h mass:  de. + D . — - + K. . ,(e. - e. J + K. . . . ( e . - e . . . * "  constant.  1-1'  M+1  i  s  i+l)  T  m  T  e  = T. - T.  (20)  . where J.  = Moment o f  i n e r t i a of the i - t h  rotational  the i n e r t i a o f a s e c t i o n o f s h a f t . D. = C o e f f i c i e n t o f speed s e l f K.j=  c o n s t a n t of the s h a f t  Units:  N-m  6. = Angle o f rotation  it  This could include  kg-m  damping.  2  Units:  N-m-s/rad  connecting masses i and j .  r o t a t i o n o f the i - t h mass, measured in the d i r e c t i o n o f in steady s t a t e  T T = The a p p l i e d mechanical rotate  Units:  (viscous)  Spring  mass.  from a f i x e d r e f e r e n c e .  Units:  rad  input torque on the i - t h mass, tending  in the + 0 d i r e c t i o n .  Units:  to  N-m  e . T. = The e l e c t r o m a g n e t i c torque on the i - t h mass, tending to r o t a t e i t in the - 0 di r e c t i o n . Equation G. and i  (20) may be expressed as  two f i r s t o r d e r d i f f e r e n t i a l equations  6. chosen as s t a t e v a r i a b l e s . i  The f i n a l  result  is:  de.  (21-a)  i  dt -I-  d6.  dt  D. i  K..  K  1  i• the masses,  model of the mechanical  T  M  E  T. - T  i  With " i " being summed over a l l mathematical  with  equations  system.  1  (21-a,b)  d e s c r i b e the  16.  3.  STEADY  Electrical  •3.1  The that  used  in  step  solution  STATE  INITIALIZATION  initialization implementing  code  the  of  used  Method  I of  t e c h n i q u e and t h a t network  is  (zero sequence  terms  and w h e r e  ref.  and b a s i c a l l y  solve  i n i t i a l i z a t i o n , only  The a s s u m p t i o n  With  EQUATIONS  f o r the e l e c t r i c a l [5], since  o f Method  b e e n c o n s i d e r e d , w h e r e no n e g a t i v e  zero)  OF MACHINE  Variables  the e x t e r n a l purposes  -  I employ  on t h e g e n e r a t o r  stator  voltages  side  lines  generator  of equations.  state  conditions  are continuously  currents  a r e f o u n d by t h e 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  solution  of  network.  These  at  currents  i  transformed When  zero  t o d - q - 0 components the machine  sequence q u a n t i t i e s  From e q u a t i o n s  is operating  are zero.  di ,  di  __i =  __s.  dt  dt  (12)  di ,  =  dt  D'D  =  the rotor angle  in steady  state,  the steady  state:  di_  n  _P_ dt  with 0  In dX  =_ i  a n d (13) R  because  present  identically [14]. the  fnom a  , 3  be  are  t=0,  Program  For have  is  transposed time  stator  phasor  i, , and i D  is  of  are present.  transformer  the transformer  s p e c i f i e d by t h e u s e r  the e n t i r e e l e c t r i c  a Thevenin equivalent  steady  of  was  t h e new t i m e -  and z e r o s e q u e n c e c o m p o n e n t s  c o r r e c t where a Delta-Wye  the transmission  both  t h e same s y s t e m  balanced  equations  cannot C  still  damper  unknown.  c u r r e n t s and  di  Q,= _g= 0 dt  (22)  dt  '  v  (22):  V Q  =  0  Vg  (23)  = 0  or i  D  =  i  0  Q  = 0  i  g  = 0  d " "Vd " "Vq v = -R i + a i L . i . +0x^/1 q sq udd v 2  (2k)  V  f  v  f  = -R i f  f  (25)  f  r  (26) (27)  17.  A phasor diagram may be c o n s t r u c t e d of  [15] from the known terminal  conditions  the machine  (V , I , and the angle between them), as shown in F i g . 5. 3 3 The r e s u l t a n t v e c t o r o f V + R I + jX I l i e s on the quadrature a x i s as a s a q a demonstrated  in r e f . [5]. Angle 5, the p o s i t i o n o f the quadrature a x i s ,  is  found from: lm(V 6  = Tan Re(V  If  the r e f e r e n c e a x i s  6 is the load a n g l e .  corresponds Using  a  + R I + jX I ) s a g a  (28)  + R I + jX I ) a s a q a J  to phase " a " o f the i n f i n i t e system,  the convention that the quadrature a x i s  DIRECT AXIS  then  lags the  QUADRATURE AXIS  REFERENCE AXIS  Fig.  5- Phasor diagram o f the synchronous machine in steady s t a t e .  18.  direct  axis  by 9 0 ° , t h e p o s i t i o n  of  the rotor  \  e = -6 + For  balanced  quantities  steady-state  operation,  (29)  d-q-0 coordinates  are related  to  phasor  by: i v  /3Te"  + ji , = d  q  where  is:  j6  + j v , = /3\/e~ ' J  a  q  I a n d V a r e RMS p o s i t i v e  sequence  (30-a) (30-b)  6  phasors  Referring  [5].  to the  phasor  diagram: i  d  i  = / 3 | I | s i n (y = /31 I I c o s ( y  -.6) "  (31-a)  <S)  (31-b)  = / J | v | s i n ( a - S)  (32-a)  q  v v  From e q u a t i o n  d  q  ( 2 6 ) , the f i e l d  ,  -  /JI VI cos(a  =  current, v  i  (32-b)  , c a n be f o u n d .  f  + R i s q  q  = -a  6)  X ,i , d d  (33)  /J This  3.2  completes  Mechanical In  and a l l  System  state,  all  t h e masses w i l l  rotate  at  (20)  Equation  .(0. i  Upon a d d i n g number o f  -  0.  J  i-1  the independent e l e c t r i c a l  variables  Variables  the steady  (I7"b).  K..  the i n i t i a l i z a t i o n of a l l  second t h e same  may be r e w r i t t e n  + K, . . , ( 9 . ii+1 i  algebraically  derivative  6. ) i+l x 1  of 0 w i l l  be z e r o  frequency, u i , given  by  m  in the form:  =  TTi  the n equations  lumped m a s s e s on t h e s h a f t ) ,  angular  terms  all  - T? i  expressed terms  in  D.oj  m  i  ,i = l,...,n { Q  nt»c>  by (34), K.  (where  disappear,  / ,  (34) n is  M  the  leaving:  6 = E  T  i=1  gives  balance  the  total  the opposing  The e l e c t r i c a l determined  input  total  centage  e T. i  torque,  by e q u a t i o n s  vidual  mechanical  of  ( T + D.u>)  (36)  e  1  the  turbine  total  physical  position  1  equation  (3*0  for  linear  algebraic  pletes  the  all  e  i=1  '  is  torque  mechanical be  equations  torque and  m  is  for  torques the  the in  turbines  to  the steady  r o t o r , where T  nonzero only  calculated  state.  e rotor  is  from  for  (36).  d e v e l o p e d by e a c h  the  turbines.  Knowing  turbine,  the  the  per-  indi-  assigned. the e l e c t r i c a n g l e o f  a n g l e may be  is  d e v e l o p e d by  damping  , are  torque  t h e masses  initialization  T  determined  this  1  (18).  torques,  can  m  '  nonzero only '  and  (13)  already of  = S  mechanical  ,  input  torques  Having  m  electromagnetic  The m e c h a n i c a l The  (35)  1  gives:  I This  - D.u) ) m  m  1=1  Rearranging  (T - T  except  easily  procedure.  found the  solved  from e q u a t i o n rotor, for  the  the  the  rotor,  (I7~a).  resulting  remaining  6..  the Writing  set  of  This  (n-1) com-  20.  4. THE TIME STEP - NUMERICAL SOLUTION OF THE MACHINE EQUATIONS  4.1 Choice o f  I n t e g r a t i o n 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 p r e s e n t i n g the complete synchronous with the r e s t of integration all rule  machine model, are s o l v e d  the e x t e r n a l e l e c t r i c network.  implicit  for solving  integration  f o r a s e t of s t a b l e  is  the machine e q u a t i o n s .  numerically stable  d i f f e r e n t i a l equations  (2) T h i s second o r d e r method is  simultaneously  s t a t e equations  i n t e g r a t i o n was adopted f o r the f o l l o w i n g  (1) T r a p e z o i d a l  re-  A l a r g e number o f numerical  schemes have been d e v i s e d f o r s o l v i n g  o f them are s u i t a b l e  ( 1 2 ) and ( 2 1 ) ,  but not  Trapezoidal-  reasons: f o r any s t e p s i z e ,  At,  [16].  self-starting,  r e q u i r i n g no past  history  terms except from the immediately p r e c e e d i n g time s t e p . (3) The c h o i c e o f At the network;  it  is is  not a f f e c t e d by the s m a l l e s t  time constant  in  r e s t r i c t e d o n l y by d e s i r e d accuracy of the s o l u -  tion. (4) There seems  to be no d e f i n i t e evidence that h i g h e r o r d e r  i n t e g r a t i o n methods  lead to numerical  ( 5 ) Use o f the t r a p e z o i d a l with  rule of  integration  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  s t e p s i z e and the same i n t e g r a t i o n (6) T h i s method is  simple  i n t e g r a t i o n methods  4.2 S o l u t i o n  superiority  implicit  [6].  simplifies  interfacing  Program, which uses the same  technique.  to program and is  the s i m p l e s t  of a l l  implicit  [18].  of the Machine Equations  by Newton's  For the system of s t a t e v a r i a b l e x = f(x,t)  Method  equations: (37)  21.  application  of  the  x(t)  Using  this  system  r  *f  =  *D  =  where  At ~T  equations  "  {  (  y  ^ - D {  R  (  differential  f {  !  D  1  i  "  V  +  +  r  +  i  g  D  e l e c t r i c angle  the for  rotor the  2  R  )  i  approximation:  equations  for  (38)  the  electrical  60  "  mass  0  -  "  )  T  *d v q  V  h^d^c?  "  R  -  (i  s  ~ 0 "  R  s  (  i  q  0  ~  +i  +  i  q  0  )  ("V^q* (wip +u$ )  ) +  r  q  T  q  }  ^ > >  + f  (39)  f  }  ) } g' and co t h e  mechanical  i-th  the  equations:  d' q' O'  i  v  . + At_ =j  angle  angular and  frequency  speed  by  of  the  (I7 a,b).  network,  The  _  mechanical  become:  e.)  (e. +  (^o)  " i I J - . - i i ii-I i-! " ii"l ii 1 i i i 1 i - 1 ? - "1 - V i ii-i^-i - i i - i " W V ii iVi 9  i  O  V  + —Y (-R g ( i 2  9  the  to  "  {  +  +  0is  related  the  in  + f(x(t-At),t-At)}}  x  {'Vd'VV  +  *f  = If  9  {f( (t),t)  v  ~ ^0  =  results  + ~ ' {-v ( i ,, i , i - , 0 ) 2 q d' q' 0'  y  ^0  rule  + £1  algebraic  +  = $ q  q  *  the  *d  =  \1>  = x(t-At)  approximation,  become  ^d  trapezoidal  +  At {  D  6  + K  (K  e  + K  )6  +  K  6  +  i  + T  +  K  (K  +  +  K  +  +  The on  T  m  subscript, the  -T }  (41)  6  "iV,  runs  turbine-generator  time  t-At,  have  have  been w r i t t e n w i t h  from  1 to  shaft.  For  been w r i t t e n w i t h  Equations  n,  no a r g u m e n t  (39),(40) , and  where  n  is  convenience,  " ~ " over  them,  the  number  all those  of  variables at  lumped at  t h e new  the time  masses old t  "(t)". (41)  are  nonlinear  algebraic  equations  22.  and an i t e r a t i v e s o l u t i o n in such  flow fast  cases because  x^,...,x  (1)  To i 1 1 u s t r a t e is  required. very  [19] » and, (2)  studies [20].  of  is  Very o f t e n Newton's  good e x p e r i e n c e  i t s convergence  l'--- m  ( x  x  (1)  method  T  proceeds as  (2) Using x ^ ,  power  equations:  (42)  - 0  m  follows  Set i t e r a t i o n count, the s o l u t i o n  network  = 0  )  g (x ,...x ) 2  in  method, assume a s o l u t i o n  r e q u i r e d f o r the s e t o f "m" n o n l i n e a r 1  it  is q u a d r a t i c and t h e r e f o r e  the use o f Newton's  g  Newton's  with  method is employed  i , to 0 and make an i n i t i a l  guess, x  (i)  , to  vector, evaluate g { x ^ }  and the Jacobian m a t r i x [ J ^ ] ,  d e f i n e d as:  3x,  3x  3x,  m \43)  [J] =  (3) Solve  9g  3g  3g,  3x, 1  3x, 2  3x m  the s e t o f l i n e a r equations [J  ]  ( i )  A^=  (4) A b e t t e r approximation x (5)  If  any one o f  (  |Ax^|,  i  +  -g(x  f o r the v e c t o r Ax: ( i )  )  = x  (  i  )  (44)  )  to the s o l u t i o n 1  x=x^  is then found from:  + Ax  . . . , | A | are l a r g e r than a predetermined s,  then increment the i t e r a t i o n count by 1 and r e t u r n to s t e p otherwise  (45)  the method has converged.  (2);  23-  Newton's  method is  a p p l i e d to the machine equations  the c u r r e n t s and mechanical angles at each new time s t e p . starting  point,  only one being size, At, extraneous  4.3  is  the d e s i r e d s o l u t i o n .  solutions  d e f i n e d as  will  by Newton's  method r e q u i r e s  r e s p e c t i v e l y . Taken t o g e t h e r ,  form o f equation  The e l e c t r i c a l  (4l) and r e a r r a n g i n g  contain  i n i t e r m s of  machine  constraint  equations be  6  equations,  hand s i d e o f equation (44).  (40) g i v e s :  terms g i v e s one  algebraic  For the i - t h mass:  ~~T—) ® • " o~T~~ ^ K. . , 6 . Jj i 2J; - - i i - l i - l  equations  can be expressed (17-b)  +  circumstances  Let g™ and g  (42).  and e l e c t r i c a l  they form the r i g h t  c o n s t r a i n t e q u a t i o n f o r each mass. 9? = (-TT i At  the s t e p -  method.  the synchronous  the mechanical system, e q u a t i o n  (46) i n t o  solutions,  Matrix  the v e c t o r s o f mechanical  Substituting  simulations,  not be generated by Newton's  in the c o n s t r a i n t  Beginning w i t h  transients  l e s s than one m i l l i s e c o n d and under these  Form o f the J a c o b i a n Solution  In  for  Depending on the  the procedure may converge on a number o f d i f f e r e n t  normally  to be w r i t t e n  to s o l v e  . - (K. . + K . . . ) 6 . + K . ^ , 0 . ^ , u-1 i i+l i i i +l i +l ]  the network f r e q u e n c y , u>.  the r o t o r angular  J  1  This  frequency  frequency from (40) with  as:  tvr-r  u = n 8 = n p r p At L  where the s u b s c r i p t  ' V  9 v  -  ( § + r At  r'  0 )']  r e f e r s to the r o t o r and n  v  is  the number o f  (48) '  pole  2k.  pairs  in the synchronous  electrical  . A  machine.  c o n s t r a i n t equations  - *  d  ( i  d.  V )  - ^  D  Using e q u a t i o n s  a r e d e f i n e d by:  -T-  +  t v  d  ( ,  d'  *Vq q V g ^ " - " ^ (,  4  - W  - p*d n  " *q  (  ,  d' P D ,  93 " • o o - * o ( ,  )  +  95 -  d  y  7  D  q  v  (  e  r" r 5  _  0.,  i  g  G  +  for  is advantageous  v  (  i  d' q' o'  6  *d  1  i  q  -T-"  [  i  : { R  ^  t o arrange  0  0  R  "T  g  q  0  P< D 'D I  +  V W .g  2  g  (i  ) }  }  +  r  )  +  V  s  R  (  d  i  +  y  }  (49-a)  R  s  (  i  =  q  +  ]  =  i  q  )  }  (49-b)  0  (49-c)  - 0  0  )}  (49-d)  0  (49-e)  °  =  =  (49-f)  °  (49-g)  g  method t o the a l g e b r a i z e d machine e q u a t i o n s ,  them i n the f o l l o w i n g  the r o t o r , Group 1 o f the v e c t o r " g " c o n t a i n s  tions  V  +  d ' f ' D>-*d  ( 1  s  -f-<VV f f 'f  f  +  {  +  )  +  Before a p p l y i n g Newton's it  )  e  o(i »i .i ) v ^ (i ^ )>  " * +  W ' P ' D ^ D  V  ,  d  f  96 - V v W  )  {  g^ = * ( i , i , i ) f  i  ,  q  =  ^ f "  +  ' o' r  ,  "^qVVV'V °  9  W  ( 1 2 ) and ( 4 8 ) , the  groups  (Fig. 6).  Except  a l l the mechanical e q u a -  b e g i n n i n g w i t h the H.P. t u r b i n e and ending w i t h the e x c i t e r and the  feed water pump.  Group 2 c o n t a i n s  with the d-q-<j> v o l t a g e rotor c o i l s  equations  a l l the e l e c t r i c a l e q u a t i o n s  f o l l o w e d by the v o l t a g e  equations  a n d , f i n a l l y , t h e mechanical e q u a t i o n f o r t h e r o t o r .  elements o f the v e c t o r o f independent v a r i a b l e s ,  starting f o r the  The  " x " , a r e arranged as shown.  25.  INDEPENDENT VARIABLES  CONSTRAINT VECTOR m 9l  '1  \ GROUP  m r-1  3  1  m r+1 m 'r+2 9  =  e 9l  r-1  3  r+l  3  r+2  i  GROUP 2  e  Fig.  6. Grouping of equations and independent " r " indicates rotor quantities.  The  Jacobian elements c o r r e s p o n d i n g to Group 1 a r e :  30.  — = . i-l  99? i  - w r I  3q  — • K 2J. Ii—l  (50) ^ '  I  .  ,  ~  ar  D,  i.  +  IT.-^M-I^IW!  +  l  ( 5 , )  l  m  (i + l Partial  variables,  i  d e r i v a t i v e s w i t h r e s p e c t to a l l  other  independent v a r i a b l e s  are z e r o .  26.  Taking p a r t i a l  derivatives,  the J a c o b i a n elements f o r Group 2 a r e as  -  S  follows:  L + d '  3g*  ai  3i  At 2  d 3i  9 v  e -(  5^t  r  +  n p  e )].L r  (5*)  q  t le7  t 1F9 g  =  n P  "  2  A  t  3  2  d  -  = n  —  [ P  " " " 3i V q  . + -^[R 2 s  L q  [n  _  At  ai  9  g  3i  f  992  9  ^ T  ai  v  & w  2  V  + n i )]. p r'  ( 5 8 )  L,  (59)  1  (60)  ]  q  (61)  o  2 _  =  p •  v  2  m  "f  2 .V 2  -yiM  (63)  n  q  g  r At 2 L  ^^[J_,  w n  92  ai  3i_ Q  r  -(  6_  p r  +  e  " 2 -  _ 1 3G  d  1e—  %  2 o e  9 v  3v +  q  2  3i  At  d  39o — ^  ai  (56)  f  g  9g  ai  M  \  D  ag* ai  17  V  9  1  1  9 i  8i  (55)  3i  e 9g  f  9 g  lp  (53)  3v  At 2  3g| '  8  ]  +  q  o  3i  IR  t  q  e .39],  9  A  2  L  d  9 i  _ q _ _ 3© d V  j J  ^  ( 6 i | )  93  9  393 3i  At 2  9 v  0  (66)  q 9 v  +  1;  ^  derivatives  in Appendix  3i  +  of v  d >  n  J T  (  ]  v , v  Q  with  respect t o  i , d  i , q  i  Q  are  6  7  )  cal-  I).  U  *lR'  +  (68)  2 "f  "  f  9c  3  3i  5  D  -  U  D  ^ 2 R „D  +  5  3i  3 i  d  994  3g|  3i  3 i  D  .  H  d  3  9  3i  '  9  -  L  g  • ^  39y  3i  3i q  q  96  ai  rotor:  d  (69)  =VZJl fM , M  = f  (70)  (71)  f r  . At  39^  9  Ror the  (65)  o  3i  o  (The p a r t i a l culated  3V  q  >«;' 3 i  At 2  9  9  3 i  R  2  v  (73)  g  2  q  6  — = n  M  (75)  28.  At jI V r  K - V  l  + K  (e  rr+1  V  +  r  §  K  r r + 1  +  U  {  D  r r  (e  r + 1 +  q  '  rr-l  +  (  T  torque,  (  e  K  r-l  6  +  9  r-l  )  +  m , ;ra )-(T; ;) ( T ^ T p  r + 1  Electromagnetic  Sp"  -  rr- 1  ( K  r  e  r + 1 +  V  r  e  r + 1  +  ) -  r r 1 >' <  e  +  + 6 r  r> (76)  =o  }  , i s known from ( 1 8 )  r>  §  +  + K  "  K  (T^  rr-1 < r-l^r-l > e  ^)  +  <  +  +..n [* (i ,i p  d  d  1 !  K  rr-1  ,i )-i D  q  +  K  rr  +  1>*W  + * -i ] d  q  (77) - W  '  Substituting  ^  '  W  d  for 6  -  +  V ' d  +  ]  °  = }  from (41) g i v e s :  r  - rr-l r-l -  17/  K  6  ( r -1  +  K  rr l>' r " " W r + l  + K  9  r  +  2 +  " *q  W d ^ f ' ^ ' q  f j + n p  +  ^ r r - l V l  {  q  d  -  i|i  +  ( K  q» Q' g ' d  ( ,  ,  ,  rr-l rr l> ' K = 0 + K  +  • ! . ] }  q d  )  ,  "  ]  "  }  K  [  A t  rr l r l i  +  +  D  n  +  T .  "  ^  "  ^  1  ^ " " r  (78)  Rrom ( 7 8 ) , the l a s t row o f the Jacobian m a t r i x i s :  3 g  9  m r  41  6 _ r  1  2J  r  K , 'Vr-l  (79)  m  ae  = - ^ K 2 J rr 1 r  +  a ; 9  ,  [ L  3i  (80)  K  r+1  d  2J jp  _ ^ ]  d q  L  r  (81)  q  v  m  ^r_= 91  9  q  qd  y  r  m  ^ m  i  3i  f  Lj]  [  pd  2J r  D  2J  p V2 f q L  r  (82)  29.  3  3  9  3 g r  i 4  "  9 i  r g  At"  r All  2 J  r"P  r  other p a r t i a l s  are  r  zero.  The nonzero s t r u c t u r e of the J a c o b i a n , The reason f o r a r r a n g i n g  the equations  Jacobian becomes p a r t i t i o n e d  as  in F i g . 6 is  i n t o two d i s t i n c t  parts.  independent p a r t , d e r i v e d from the mechanical T h i s p a r t comprises  12% o f the m a t r i x  winding machine model, that  does  case.  lumped masses  mass system,  is the  is as  l a r g e as  the t i m e rotor. seven-  possible.  comprises o n l y 2 8 % of the m a t r i x  the s i z e of the time-independent p a r t of  case o f a o n e - ( i n f i n i t e )  The f i r s t  equations without  i s , when the m a t r i x  As the number of  now more e v i d e n t , the  in the case o f an eight-mass  The time-dependent block on the lower r i g h t for this  [ J ] , is shown in F i g . 7-  in the model d e c r e a s e s , [J] u n t i l ,  f o r the  so  limiting  t h e r e is no time-independent p a r t  at  al 1 . For s i m u l a t i o n  purposes,  the s t r u c t u r e of the J a c o b i a n .  is  process  the e q u a t i o n f o r the r o t o r ,  for solving  saved by e x p l o i t i n g need  i t e r a t i o n s o r from one time step to the n e x t .  reduced d u r i n g LU t r i a n g u l a r i z a t i o n of  stitution  is  The time-independent p a r t , o b v i o u s l y ,  not be updated d u r i n g s u c c e s s i v e A l s o , by w r i t i n g  some computer time  the AX..  last,  the number o f  eliminations  the m a t r i x p r i o r to the back  sub-  CD: < T tn  < (O  <  <  -ti  II  O  < II  <  Xl  II  <  Q_  <_ OO  t—v n c _  CD : CD oo ll  II  II  r -  3' .7+O O O 0--:-s (D — QJ ^ TJ -I  O  D_  C  n  3  CD —  m.  (D -i O in  lO  zy ~\ rt  C  i o  — 3  3 to  "J N  QJ in in  rt C -i (D  oP  o m -D m z o m z H  ^>  uP ~P  X X  X X X X  rt  3" in (D  —' Q) -I  in  O -t  (T> < (D  rt  I  (D  3  oj  3  3  a n — o 3 3 cr  rt  QJ  IQ  rt IQ  — X  CD  •  —  QJ  -I  3  CO  coCD  3  3 QJ  -1 r t QJ -1 rt — O X  H c—•  rn o m ~o rn z rn z H  O  o  LO  X X X X X X X XX X XX X X X X X X X X XX X X XX X X X XX  Ul  c_  <_ .  c_  —* CD : CD : CD CD : a>' UJ -f ro —* II  II  UJ  II  ll  ll  X X XXX XXX XXX XX  xl X X X X X X X  XX  1  UJ  • TIME DEPENDENT  TIME INDEPENDENT  o  31.  5. STEADY STATE AND TIME-STEP SOLUTION PROGRAMS  5.1 P r e l i m i n a r y Remarks About the Programs The f o l l o w i n g nous g e n e r a t o r model (1)  is  the e l e c t r i c a l  i n f o r m a t i o n must be s u p p l i e d by the user to be used in the T r a n s i e n t s  at the g e n e r a t o r t e r m i n a l s and the  angles o f the t h r e e phasor v o l t a g e s w i t h  72 f *3  TjM , r  (5)  resistances  the number o f poles the s h a f t ,  R , R, s f  L,.  f  Lr,,  D'  R^, R_, D Q  £  L  sections,  f o r the s h a f t  (8) damping c o n s t a n t s ,  q'v2q  ,  L  g  .  L  n>  Q  and  L  rt  in  h.  0  lumped masses on  kg-m . 2  in N-m.  in N-m-s/rad.  (9) the percentage o f the t o t a l  mechanical  torque c o n t r i b u t e d by each  stage.  In a d d i t i o n to the above, small  A small  f3  reference.  and the mass number c o r r e s p o n d i n g to the r o t o r .  (7) s p r i n g c o n s t a n t s  Transients  ,  in the g e n e r a t o r , the number o f  in  resistances  r e s p e c t to the  and R in ft. g  (6) i n e r t i a s of the lumped masses,  turbine  Program.  frequency.  (2) the peak l i n e - t o - n e u t r a l v o l t a g e  (4) the f i v e  i f a synchro-  must be s p e c i f i e d so t h a t ,  z e r o sequence and p o s i t i v e sequence  in the steady  state c a l c u l a t i o n s ,  the  Program can c o n v e r t the g e n e r a t o r to a Norton e q u i v a l e n t c i r c u i t .  but f i n i t e r e s i s t a n c e  (»0.001ft)  is  required across  each c u r r e n t source.  T h i s program c o n v e r t s the two r e s i s t a n c e s  to s e l f and mutual q u a n t i t i e s which  are subsequently  r e s i s t a n c e s o f the machine d u r i n g  steady s t a t e  i n c l u d e d with the s t a t o r  initialization.  Even the s m a l l e s t  may be l a r g e compared to the a c t u a l  coupling  armature r e s i s t a n c e .  problem, the zero and p o s i t i v e sequence r e s i s t a n c e s each o t h e r , f o r  i n s t a n c e , 0.001ft.  resistance permitted To overcome  this  should be s e t equal  to  By s p e c i f y i n g an armature r e s i s t a n c e equal  32.  to the true value minus absorbed  i n t o the  0.001ft, the f i c t i t i o u s r e s i s t a n c e disappears when  stator.  The program assumes t h a t , high p r e s s u r e  t u r b i n e appears  first,  stage and so on down the s h a f t . model  the mechanical  synchronous  machine.  corresponding  in l a b e l l i n g  It  system but r e q u i r e s  the user does not wish  the e l e c t r i c a l  large  a t w o - p o l e , one-mass  r o t o r mass having  Q u i t e o f t e n , the m a n u f a c t u r e r ' s i n c l u d e : X,, d  The f u n c t i o n a l  r e l a t i o n s h i p between these constants  eters  and mutual  1  inductances  is  method, d e s c r i b e d in r e f . '  1  i n s t e a d o f the o p e n - c i r c u i t v a l u e s .  d i r e c t axis  can be o b t a i n e d with  parameters  T ' , T " , and qO qO  partly nonlinear.  Sometimes  [5].  1  frequency.  and the d e s i r e d s e t  r e q u i r e d by the program may be e x t r a c t e d from t h i s  Newton's given  and s e l f  1  rotor  f o r a machine  rotational  T ' , X , X , X", dO' q' q q  the  system w i t h the  data o f the e l e c t r i c a l  of the machine  resistances  X ' , X'!, T ' d' d dO  constant  to  d e s c r i p t i o n of  to mass number 1, the t i m e s s t e p program s o l v e s  w i t h an a r b i t r a r i l y  the  f o l l o w e d by the i n t e r m e d i a t e p r e s s u r e may be that  By s p e c i f y i n g  the masses on the s h a f t ,  data by  R. f c  of  The paramapplying  T' T'J, T ' , and T " are dl d ' q' q  The d e s i r e d time constants  for  the  [21]: X  T  io " 7  1 T  Xd  i  <>  'd  (87)  86  K 'dO  Y 1 1 X  Extension  to the quadrature a x i s If  rature axis mode 1.  L  is  s e t equal  is  obvious.  to z e r o ,  and the program s o l v e s  d  t h e r e is  no " g " winding on the quad-  for a six-winding  electrical  machine  33.  5.2  Steady  State  - Subroutine  The f o l l o w i n g steady  state  subroutine I  steps  i n d i c a t e the b a s i c  flow o f events  i n i t i a l i z a t i o n of machine c u r r e n t s and mechanical  during  the  angles by  "MSTEA .  Initialize (1)  "MSTEA"  11  Electrical  Calculate  real  Equations and imaginary  sequence c u r r e n t s  parts  of p o s i t i v e , negative,  from real and imaginary  c u r r e n t s given by the T r a n s i e n t s (2) The magnitude of n e g a t i v e degree o f unbalance  parts  and z e r o  of the terminal  Program.  and zero sequence c u r r e n t s  in the system.  indicates  Both should be very small  the com-  pared to the p o s i t i v e sequence c u r r e n t . (3) C a l c u l a t e  the real  and imaginary  (k) F i n d the angle of  of  lm(V 1  Re(V  a  + R I + jX I ) s a q a  a  + R I + jX I ) s a q a  £ = 6 j and i t s  angle,  3 ground v o l t a g e Solve  for  voltages: d  ^ l  =  v  V  a l  5  i  n  '  " )  a  6  = /3|V  | cos (a - 6) 3 and i t s a n g l e , y .  C]  (8) Compute magnitude of  I 3  Solve  for currents: i i  (10).Calculate  d  a.  V  is 3  f o r phase " a " .  v  (9)  + R I + jX I s a q a  J  (6) Compute the magnitude o f V  (7)  V  the quadrature a x i s , <5:  & = Arctan  (5) Rotor angle  parts  = / J | 1 | s I n ( y " <5) a  = fi\ l | c o s ( a  the f i e l d c u r r e n t , i ^ .  Y  - <S)  the RMS  line-to-  34,  v  All (11)  q  +  R  i  - (iiL.i ,  s q  d d  o t h e r c u r r e n t s are z e r o .  Form P a r k ' s  inductance m a t r i x ,  f-pl*  (12) Transform v o l t a g e s and c u r r e n t s (13)  II  Initialize  Initialize  the v e c t o r o f f l u x  Mechanical  to phase c o o r d i n a t e s  linkages,  ip =  [L^ ] i .  Equations  (1) C a l c u l a t e e l e c t r o m a g n e t i c torque on the r o t o r plus speed s e l f - d a m p i n g  input torque to each t u r b i n e .  t u r b i n e torques equals  the torque c a l c u l a t e d  (3) Form the c o e f f i c i e n t m a t r i x ,  [K]  in  The sum o f  the  Il-(1).  , of the mechanical  equations,  hand s i d e v e c t o r g. = T ? - T? - D.cj , again 1  m  omitting  rotor.  (5) Solve  the s e t o f  the i n i t i a l  l i n e a r equations  mechanical  [K]6 = g f o r the 6..  These are  angles at time t = 0.  (6) C a l c u l a t e the mechanical HISTl(i)  all  the r o t o r e q u a t i o n .  (4) Form the r i g h t the  the sum of  torques.  (2) A l l o c a t e steady s t a t e  omitting  for printout.  history  vectors:  = 2 u> + (^- + ji-)8. i m  = | T t-K.. i - 1 ( -1 + K.. . J e . 2J. i-1 M-1 1 i+l 1 1 (7) Return to the T r a n s i e n t s Program.  HIST2 (i)  e  +  5.3 Time Step - Subroutine  K  •<.._,, 9. 11 + I  TT +  1  T?} 1  "GEN"  The t i m e - s t e p s o l u t i o n o f the synchronous o f the f o l l o w i n g  -  i +1  machine equations  consists  steps:  (1) P r e d i c t new mechanical  angles  8. = B. + 9.At and new c u r r e n t s  i = i.  35.  ( " „ " means  from the p r e v i o u s  these p r e d i c t i o n s (2)  Modify ing  (3)  stator  time s t e p ) .  In steady  state  operation,  the small  external  are c o r r e c t .  resistance  v e c t o r to absorb  coupl-  resistances.  Compute the time-independent constants matrix, actual  "JACOB", and the v e c t o r " G " . F o r t r a n c o d e . i n Appendix  appearing Terms  in the J a c o b i a n  in quotes  appear  in the  III.  (k) Compute the time-independent elements  (mechanical  terms)  of  the  Jacobian. (5)  Calculate  the r o t o r e l e c t r i c a l  past h i s t o r y  angle and i t s  angular  and the p r e d i c t e d new mechanical  (6) Compute the time-dependent f a c t o r s  appearing  frequency from  angle. in the J a c o b i a n and the  vector " G " . (7) C a l c u l a t e t e r m i n a l  voltage  and f l u x  linkage  vectors  to the p r e s e n t e s t i m a t e of c u r r e n t s and r o t o r (8) F i l l (9)  in the elements o f  C a l c u l a t e the v e c t o r s terms  (10)  corresponding  angle.  the time-dependent p a r t of the J a c o b i a n .  of the mechanical  part  that w i l l  become h i s t o r y  hand s i d e v e c t o r " G " of mechanical  and e l e c t r i c a l  in the next time s t e p .  Calculate  the r i g h t  equations. (11)  Solve of  f o r the c o r r e c t i o n terms, AX,  linear  in Newton's method from the  equations:  [j] .AX = - i (12)  Update the c u r r e n t s and mechanical i j  L . f - 0  eK j  k  +  A  i  .  + A0 . j J  j  angles:  set  36.  (superscript (13)  If  "k"  indicates  i t e r a t i o n number; k = 1 , 2 , . . . 1 0 )  more than ten i t e r a t i o n s have been performed, e x i t from the p r o -  gram w i t h a warning message. (14)  If  any of  epsilon (15)  |AX.| are l a r g e r  The procedure is not  converging.  than p r e d e f i n e d t o l e r a n c e v a r i a b l e ,  (=10 ) , go to (5) and i t e r a t e once more. 9  Since convergence is  successful  this  p o i n t , update terminal  voltages,  flux  (16)  Store a l l  quantities  (17)  Convert d - q - 0 c u r r e n t s to the phase domain to be i n j e c t e d i n t o the external  (18)  linkages,  at  and mechanical h i s t o r y  vectors.  needed f o r the next time step  solution.  network.  Return to the T r a n s i e n t s  Program.  5.4 Flow Diagram of S o l u t i o n by Newton's  Method  Under the c o n t r o l o f 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, time step program s o l v e s  the synchronous machine equations  step o f the s i m u l a t i o n .  F i g . 8 shows the flow diagram o f the s o l u t i o n  a l g o r i t h m by Newton's machine programs  the  d u r i n g each time  method, with the i n t e r f a c i n g between the t r a n s i e n t s  in phase c o o r d i n a t e s .  and  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  Fig.  8. Flow diagram o f s o l u t i o n by Newton's method.  38.  6. TESTING OF THE SIMULATION METHOD  The main purpose o f  running t e s t cases was not to v e r i f y the s y n -  chronous machine model, s i n c e t h i s  has a l r e a d y been done by o t h e r s  r a t h e r , to c o n f i r m the p a r t i c u l a r s i m u l a t i o n method chosen. done in v a r i o u s 6.1  Initial  [5], but  The t e s t i n g was  stages.  Tests  As a p r e l i m i n a r y t e s t , a very simple network was chosen so the T r a n s i e n t s Program could be r e p l a c e d by a Thevenin e q u i v a l e n t The network, shown in F i g . 9, c o n s i s t e d o f phase r e s i s t a n c e s o f nected to s i n u s o i d a l torsional  that  circuit.  10.fi c o n -  v o l t a g e sources w i t h amplitudes o f 5000. v o l t s .  Shaft  e f f e c t s were n e g l e c t e d by m o d e l l i n g the mechanical system as one  arbitrarily  large  r o t a t i n g mass.  were those used in Example k of to 10,000. v o l t s .  The e l e c t r i c a l ref.  [5] and peak terminal v o l t a g e was  The machine was s i m u l a t e d  hundred 200-psec. time s t e p s .  All  parameters o f the machine  d-q-0  set  in the steady s t a t e over f i v e  c u r r e n t s were constant with time.  However, when the network was modelled in the T r a n s i e n t s Program, very s m a l l , stable o s c i l l a t i o n s  about the steady s t a t e values were o b s e r v e d .  GENERATOR  R = 10.fi  The a m p l i -  INFINITE BUS  V\AAA VNAAA 10000.10°\I  0  F i g . 9. Simple network c o n f i g u r a t i o n .  0  5000.A0°V  39.  tildes were v e r y 0.0004  percent of  c o u l d not  be  6.2  was  achieved  gram c a u s e d  Single  and  the steady  entirely i n one  one  convergence  of  ref.  out o v e r rents, cent of  error  as  field the  hundred  value.  The cause  the machine In  the  of  program  first  reached a f t e r  it  this and  ease,  interface with  amounted  to  oscillation its  i n f l u e n c e on  the s o l u t i o n  the T r a n s i e n t s  two  iterations  at  t h e same e l e c t r i c a l  parameters  and was  before.  is  This  test  case  1ine-to-ground shown  time s t e p s  in  with  10.  a step  o b t a i n e d by M e t h o d  is  fault  Fig.  c u r r e n t and e l e c t r o m a g n e t i c results  the " d " c u r r e n t ,  each  was Pro-  time  step.  Fault  had  configuration  four  in  of  negligible.  [ 5 ] , where a s i n g l e  The n e t w o r k  state  t o be  Line-to-Ground  i n f i n i t e mass  the case  i t e r a t i o n , while  The g e n e r a t o r as  in  t r a c e d t o any  the s o l u t i o n always  small,  that  is  applied  size  of  to  Example k phase was  500 y s e c .  I of  [5].  This  slight  "a".  carried  Phase  t o r q u e were always w i t h i n  cur-.  0.02  per-  disagreement  INFINITE BUSBAR ^  TRANSMISSION  63.2mH  in  The s i m u l a t i o n  3 SINGLE PHASE UNITS 1 3 . 8 / 1 3 0 . 8 kV  L =  used  modelled  LINE  L  =  0.22H  L"  =  0.096H  13.95A-30°kV 137-23/--20°kV FAULTED BUSBAR  Fig.  10. N e t w o r k c o n f i g u r a t i o n single line-to-ground  for fault.  40.  between small with  the  test  and  discrepancy Method  cases.  II  in  of  the  the data  [5] was  Method  6.3  IEEE S u b s y n c h r o n o u s  The g e n e r a t o r  case  for  the  a c h i e v e d when of  was  probably  caused  two s i m u l a t i o n s .  Exact  identical  Method  is  simulation  Resonance  provides  II  was  test  d a t a was  available  of  electrical  perience  at  circuit,  properly  t u n e d , was  as  observed  in  configuration  is  rather  subsynchronous  t h e Mohave  in  to  by  very  agreement  used  the  a  in  author  both but  plant to of  11 and  Case test  of  a one-mass  The  resonance  found  Fig.  than  model.  generating  the a n a l y s i s  shown  Benchmark T e s t  a much more s e v e r e  now a s i x - m a s s  a seven-winding  those  simulations  l>)  This  the  used  (The c o m p i l e d v e r s i o n  not  having  reference  solution  mechanical  IEEE b e n c h m a r k  [22]  was  based  in Southern  produce  the  Nevada.  t h e same  the a c t u a l  on  system.  the mechanical  test  method.  system, case  the a c t u a l A simple  transient  of  B  RLC  the  network  generator  INFINITE BUSBAR  X .-j:v371 c  UNITY  I—  VOLTAGE  close open  Fig.  11. E l e c t r i c network f o r the of  subsynchronous  simulation  resonance.  ex-  problems  The e l e c t r i c  model  for  t=0. t=0.075s  41  HP  )  ' )) |"p  Fig.  is  shown  ref. was  in  Fig.  represented  Mechanical  12.  the Thevenin e q u i v a l e n t  then back  to  per  unit  units  time  crossed of  t =  0  zero. ysec.  500  and  Fig.  ends  the s h a f t  oscillations shaft did  failure  occur at  ment w i t h Edison  of  section  a  this  in  after  m u l t i p l i e d by  system.  point.  The  In  of  was  rotor,  in  was  in  and  at  bus  This  Program  Fig.  in each  to and  "c",  the  would  the a c t u a l  in  the  growing  shaft  very  Southern the  The  step  between  eventually  case,  phase  simulated  torque  constant.)  11  a time  d i f f e r e n c e between  and w i t h  re-  routine  with  the s h a f t  shaft  in  network  B in  current  phase  d e v e l o p e d by  and E l e c t r i c Co. [22],  given  found.  presented here are  the program  the  system.  seconds  0.4  the s p r i n g  fact,  is  time s t e p  the  the angle  of  while  unit  the  as  current  section  results  per  applied  sec.  simulation  last  units  data  from the T r a n s i e n t s  f a u l t was  0.075  system.  machine  the s o l u t i o n  the s i m u l a t e d  real  Gas  the b e g i n n i n g  torque equals  o b t a i n e d by  and P a c i f i c  at  the  the  the  network  the g e n e r a t o r  (Shaft  torque  in  those  shows  in  the  removed a f t e r  t o r q u e on  and e x c i t e r . of  13  physical  Program  three-phase  The p e r i o d o f  electromagnetic rotor  then  in  ^)) EXC  ^)\QT0R\  the s h a f t  of  quantities  A simultaneous  of  \  and m e c h a n i c a l  simulated  the T r a n s i e n t s  physical  J\LPB  model  The e l e c t r i c a l  12.  in  be c o n v e r t e d t o  at  f)\ LPA ^  T h e m a c h i n e was  [22].  quired  \  lead  to  damage  good  agree-  California  results  given  in  [5].  hi. V  i  r  0.16  0.24  0.2  TIME (SEC) Fig.  13.  Subsynchronous  resonance,  electromagnetic  torque  EXC-rotor  torque  shaft  on  phase the  c current  rotor  (lower).  (upper),  (center),  and  43.  EXAMPLES OF PROGRAM A P P L I C A T I O N S  7.  Examples cations same  of  the  machine  resulting  synchronous  is  used  from  Sustained  in  the  phase  "a"  taneous  and  low p r e s s u r e was  Faulty  terminal  Fig.  0.5  voltage and  point  26.0Z0  ground the  torque  first  example,  120° e l e c t r i c a l  the  Fig.  and  a step  or,  and  two  to  torques  and p a r t  of  chapter.  simulate  generator.  the  high  section  of  simul-  time  when a  t =  Delta/  torque  span  the  of  0,  simul-  the  between  The t i m e  a  At  side  on last  the  sim-  psec.  on e i t h e r at  voltages  faulty "a"  angle  cases  excitation  For normal  angle  phase  error  both  used  no-load  phasors  these  shaft  electromagnetic  300  The  Generator  having  generator  In  of  appli-  Generator  rotor.  size  this  1ine-to-1ine)  the network.  when  in  14 was  (RMS  Unloaded  section,  9 = 240° e l .  used  the s h a f t  the  of  The g e n e r a t o r  unloaded  kV  0  in magnitude  (el.),  a comparison  Unloaded  in  voltage  this  Program.  simulated  generator  In  the T r a n s i e n t s  o c c u r r e d on  the  In  is  in  to  or  angle  an  a generator  occurs  possible  is  initially  an  the  permit  study  connected  nitude  error  of  in  some o f  disturbances.  is  synchronization  by  an  (LPB)  with  case,  Faulty  "a"  to  stage  correspond  angle.  to  for  15 shows  seconds  this  model  shown  was  Synchronization  the network  the  on  the mechanical turbine  In  nection  fault  fault  rotor  of  Fault  diagram  three-phase  the  7.2  types  voltage  transformers.  illustrate  the examples  terminal  Wye  ulation  machine  Three-Phase  three-phase  to  IEEE b e n c h m a r k  The s y s t e m taneous  presented  from d i f f e r e n t  the network  7.1  are  the  and  synchronization,  side  of  time o f  do n o t  3 =  leads  120°.  the magnitudes  the  In of  con-  switching.  correspond  synchronizations voltage  rated  are  in  simulated.  the  system  the  second,  the  mag-  phase the  voltages  kk.  A4 CLOSE  AT  t=0  CLOSE  AT  t=0  CLOSE 2 6 . Z_0°kV  Fig.  o.o  o.os Fig.  14.  o.i 15.  26./538.7 SINGLE PHASE  System diagram  o.is  0.2  Electromagnetic LPB-rotor  kV UNITS  for  TIME  three-phase  0.25  (SEC)  t o r q u e on  mechanical  AT  t=0  0.3  0.35  rotor  torque  fault.  (upper),  (lower).  0.4  and  0.'45  0.50  h5.  at  the  connection  voltage would  at  the  flow  connected  the  infinite on  seconds  system,  bus  was  quite  120° for  el.  quite well  from  Figs.  17  The el.,  produce  benchmark  in  having  time  Fig.  mechanical  were  The  16. per  adjusting  used  no  the  current  previously  transmission  unit  torque  and  line  to  per  X^=0.50  in  in  shows  17 the  Fig.  over the  a  shows  span  simulated  LPB-rotor  18  time  unit  shaft  the  of  electro-  for  an  error  corresponding  el.  (motor  operation).  ref.  [23].  Two c o n c l u s i o n s  may  error  1 2 0 ° and  240°  given  simulated  Fig.  operation).  of  by  synchronization,  The g e n e r a t o r  XQ=1.56  steps.  synchronizations  different  torques  synchronization  are  different  MVA w i t h  proper  network. shown  achieved  These  plots be  com-  drawn  with  and e l e c t r i c a l  with -  both  error  torques  angles  mechanical  angle  than  120°  faulty  of  240°  and e l e c t r i c a l .  produces  more  severe  synchronization  of  el.  These  780  with  was  kV).  538.72  angle  faulty  A faulty  240°  as  This  18.  mechanical  somewhat  the  those  and two  and  (generator  to  that,  long,  ysec.  an e r r o r  pare  (2)  and  same.  synchronizations  witha300  angle  (1)  generator  faulty  the  so  the  torque  torques  bus  to  magnetic of  were  8 9 2 . 4 MVA a n d The  0.5  infinite  between  was  (based  points  five  essentially machine  masses  generator  with  and  the  same  network.  on  the  an  8 9 2 . 4 MVA  conclusions The  Brown  turbine-generator rating  and  reached  Boveri shaft  having  in  [23]  generator  compared six  to  rotating  for  a  was  rated  the  IEEE  masses.  46.  26.^0° kV 3, SINGLE  TRANSMISSION  PHASE  LINE  Z =162.6 + J 5 0 7 . 3 P.  UNITS  Q  26/538.7 kV  Z^.5  + j 162.6 n  6 6 6 380.95^30 -6 kV 0  (1) 6=120° (2) 6=240° Fig.  16.  System diagram  for  faulty  synchronization.  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  This models  of  the  thesis  has  described  synchronous  magnetic  Transients  complete  set  RECOMMENDATIONS FOR FURTHER  machine,  Program,  the e l e c t r i c a l as  i t was  machine e q u a t i o n s  simultaneously  The s y n c h r o n o u s  m a c h i n e model was  interfaced with  of  terminals. studies solution  a Thevenin equivalent The  where  r e l i a b i l i t y of  c o m p a r i s o n was  techniques. of  the machine  was  to  draw  that  possible  agreed with There  Some o f  these  (2)  improvement o f variables  of  program.  made by a  a  the U.B.C.  solving  the T r a n s i e n t s  method has of  the  network.  Program  from the  machine  been v e r i f i e d  other e n t i r e l y  i n t e n d e d as faulty  for  Electro-  the e l e c t r i c  the network  results  were In  about  examples  of  synchronization  in  test  different possible  study  t h e e l e c t r o m a g n e t i c and s h a f t  it  torques  manufacturer.  possible  to  saturation  areas  of  additional  investigations.  of  the  steady  in  of  i n i t i a l i z a t i o n of  the  electrical  and z e r o s e q u e n c e c o m p o n e n t s  on e l e c t r o m a g n e t i c and s h a f t  the s h a f t  by m a n u f a c t u r e r s , (4) e x a m i n a t i o n  state  in  the  [14],  the e f f e c t s  damping  effects,  include negative  initialization  ical  the  with  torsional  are:  inclusion  study  to  studies  a number o f  (1)  (3)  made  conclusions  studies are  the s o l u t i o n  Two t e s t  applications  c i r c u i t of  to  technique  of  by means  and s h a f t  added  and a new s o l u t i o n  RESEARCH  system,  which  is  not  torque of  yet  being  mechan-  measured  and,  the s h a f t  torques  during  c o r r e c t and  should  prove  faulty  ^synchro-  nization. The s y n c h r o n o u s studying  shaft  torsional  generator  model  and o t h e r e f f e c t s  in  t o be v e r y  synchronous  machines  useful  in  following  transient second.  disturbances,  where  the  time  period of  study  may  be as  long  as  one  50.  REFERENCES  1.  V.  Brandwajn  and H.W.  Dommel,  Electromagnetic Transients  2.  IEEE  PES  IEEE  Committee  vol.  K.  Syst.,  Meyer,  E.H.  V.  tromagnetic  B.C.,  W.S.  Meyer,  7.  E.W.  Kimbark,  Dover, B.  and  Machines:  R.G.  Memorandum,  Paper  V-D,  !the  subsynchronous IEEE  December 23,  Ph.D.  thesis,  "MANTRAP,  presented at  J u n e 2-4,  Program  of  Trans.  1977Machine  t h e 9th  and  PICA  1975. for  the S i m u l a t i o n  University  to  Memorandum,  Stability:  1956  Harley,  Application  of  and power s y s t e m s " ,  G e n e r a t o r Models  System  (reprint  at  of  British  of  Elec-  Columbia,  1977-  Transients  1968  Adkins  April,  "Power  an  19-23, 1976.  the study  J . J . LaForest,  La.,  "Synchronous  Transients",  6.  and  Program",  New O r l e a n s ,  Brandwajn,  Vancouver  Program  Lenfest,  Transients  76359-0 . . p r e s e n t e d  July  for  With  PAS-95, p p . 216-218, J a n . / F e b . 1976.  Transients  Carlson,  Ore.,  machines  G e n e r a t o r Models  Paper A  bibliography  Power A p p .  Conference,  8.  "A  rotating  Network  5.  Report,  Program",  Portland,  resonance between  3. W.S. 4.  Summer M e e t i n g ,  "Interfacing  March  10,  1974.  Synchronous  Machines",  Theory of  Alternating  New  York:  ed.).  "The General  Practical  Problems",  Chapman  Current  and H a l l ,  London,  1975. 9.  G.  Gross  and M . C .  Simulation published 10.  in  "Synchronous  the Computation the  IEEE T r a n s ,  on  Brandwajn  of  Electromagnetic Transients",  R.H.  Dommel,  Conference, Oct.  Park,  eralized  and H.W.  of  V.  arid Power 11.  in  Hall,  "Two  Methods  Reaction of  Machine  Dynamics  Electromagnetic Transients", Power A p p .  "Synchronous Proceedings  20-22, 1976,  Theory of  Analysis",  and T o r s i o n a l  AlEE  pp.  and  the  Parameters  in  Analysis  Canadian  Communications  393"395.  Synchronous Trans.,  be  Syst.  Machine of  to  Machines  vol.  48,  -  Part  July  1929,  Trans.  Amer.  I,  Gen-  pp. 716-730. 12.  Y.  Yu,  Elect.  "The Torque Tensor of Engrs.  Ill,  vol.  the General  Machine",  81, p p . 623-629, F e b . 1963-  Inst.  51.  13-  P.J.  Fenwick, "Dynamic  University 14. V.  of  Brandwajn,  tialization Transients Society 15.  P.M.  Anderson  H.W.  17.  18.  G.  and A . R .  for  Trans.  Power A p p .  Gear,  19.  20.  H.W.  21.  Syst.,  vol.  Bergen,  Syst.,  23.  Canay,  ances",  Electromagnetic Engineering  of  PAS-91, Value  Boveri  and  Stability",  1976. Stability  July/Aug.  Solutions", pp.  1972,  New M u l t i s t e p  Power S y s t e m  Problems  pp.  1977,  N.  1643-1650. Algo-  Response",  in O r d i n a r y  Cliffs,  IEEE  Integration  Dynamical  Jan./Feb.  "Numerical  J . ,  IEEE  293"306. Differential  1974.  University  of  British  Methods  with  Fortran  IV  Case  1972.  Kingsley,  and A .  New Y o r k ,  Resonance  vol.  "Stresses  Brown  an  \n\t  1974.  the Computer S i m u l a t i o n Syst.,  Transient  Englewood  B.C.,  Toronto,  IEEE S u b s y n c h r o n o u s  M.  of  McCracken,  F i t z g e r a l d , C.  Power A p p .  Machine  IEEE Power  Control  on Power S y s t e m A n a l y s i s " ,  and D.D. Wiley,  Iowa,  Class  vol.  Prentice-Hall,  W.S.  for  Within  1978  System  PAS-91,  Initial  3rd e d . , McGraw-Hill, 22.  "A  Computation  Vancouver,  A.E.  the  Units",  1976.  Synchronous  Conditions  Ames,  "Fast  Columbia,  Studies",  Press,  Sato,  Dommel, " N o t e s  Dorn  B.C.,  Dommel,"  Fouad, "Power  "Numerical  Equations",  and H.W.  t o be p r e s e n t e d a t  and A . A .  rithms  C.W.  Turbo-Generator  Meeting.  and N.  the  of  Vancouver,  U n b a l a n c e d Network  Power A p p .  Gross  Meyer,  University  Dommel  Trans.  W.S.  Analysis  Columbia,  Program",  Iowa S t a t e 16.  British  for  Summer  Response  PAS-96,  Force  Report,  Subsynchronous Sept./Oct.  in Turbogenerator Sets Review,  "Electric  Machinery",  1971.  Task  of  Kusko,  vol.  62,  "First  Benchmark  Resonance",  1977, due  September  PP-  IEEE  Trans.  1565"1570.  to E l e c t r i c a l 1975,  Model  pp.  Disturb-  435 443. _  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) i n phase coordinates. d-q-0  coordinates, equations  To begin the transformation of (19) to  (3) to (7) must be used.  Equation  (19)  becomes:  v  d q  s 5  °=  1  J l i J  i s  d s  ^°  1 + L i J  [P] y "No  (1-1)  v. d v q  with v  [PHZ/^HP]" \$  dqo s  V  (1-2)  o [Z^^]  No  V  =  ^  = [VZA;VZB|VZC]  ( 1  •  _  3 )  (1-4)  VZA, VZB, VZC and E l are column vectors of dimension three.  For  convenience,  they have been given the same alphanumeric names as the corresponding vectors in the Fortran code. 2  The following trigonometric i d e n t i t i e s w i l l be useful. 1 1  cos 0 = -j + Y cos20  (1-5)  sin 0 - \ - \  ( 1  2  cos20  cos0 sinG = y sin20 2  _  6 )  ^  2  cos 0 - s i n 0 = cos20 cos0 cos(0 - 120) = -  (1-8) (1 + cos20) + - | sin20  (1-9)  53.  cosG cos(0 + 120) = - f- (1 + cos26)  4  - -™ sin20 4  (1-10)  cos0 sin(0 - 120) =  (1 + cos20) - 7 sin20 4 /3 1 cos0 sin(0 + 120) = (1 + cos2©) - — sin2 0  (1-11)  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)  4  T  4  4-  1  i/3  1 0  .  (1-12)  sine cos(0 - 120) = 7- (1 - cos2G) - 7 in20 4 4 -/J 1 sin0 cos(0 + 120) = -j (1 - cos20) - j sin20  (1-16)  cos (0 - 120) =4-7- cos20 - -~ sin20 2 4 4  , (1-17)  s  2  1  . (1-15) 1c  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 4 4  (1-21)  cos(0 - 120) cos (0 + 120) = - i+'i--cos2 0 4 2  (1-22)  sin(0 - 120) sin (0 + 120) = - 7 - \ cos20 4 2  (1-23)  cos (0 + 120) = j - j cos20 + ~  (1-24)  2  sin  2  (0 - 120) = j + ~ cos2 0 +  sin20 sin20  sin (0 + 120) = j + i cos2 0 - - | sin20  (1-25) (1-26)  Where 0 i s the angle between phase "a" and the direct axis, measured in the direction of rotation.  54.  Define [Z 1 by  following expression:  t n e  (1-27) or  2 [ Z  P  3  !  cosG  cos(0 - 120)  cos(0 +  120)  VZA(l)  VZB(l)  VZC(l)  sinG  sin(0 - 120)  sin(0 + 120)  VZA(2)  VZB(2)  VZC (2)  VZA(3)  VZB(3)  VZC(3)  1 /2  1  1  /2  S2  cos 0  1  sin 0  /2  cos(0 - 120) sin(0 - 120)  1 /2  cos(0 + 120) sin(0 + 120)  1  /2  (1-28) The following can be obtained from (1-28) with equations (1-5) to (1-26) Z (1,1) = -| cos 0 VZA(l) + \ cos0 cos(0 - 120) VZA(2) 2  + -| cos0 cos(0 + 120) VZA(3) + j cos0 cos(0 - 120) VZB(l) + | cos (0 - 120) VZB(2) + | cos(0 + 120)cos(0 - 120) VZB(3) 2  + | cos0 cos(0 + 120)VZC(1) + -| cos(0 - 120) cos(0 + 120)VZC(2) + -| cos (0 + 120) VZC(3) 2  =  (  3  +  \  cos20  + (_ A _ 1  c  o  ) s  VZA  2  0  (D + (~\_ fl  6  \  cos20  +  sin20)VZA(2)  sin20)VZA(3) + (~ \ - \ cos20 + ^  sin20)VZB(l)  55.  + (T- T +  c o s 2 0  - ^ ? sin29)VZB(2) + (- T + -77 cos20)VZB(3)  \ ~ \ os29 - i-jr  sifL20)VZC(l) + (- -jr + -| cos2G) VZC(2)  c  + (| - i  cos2e +  sin29)VZC(3)  (1-29)  Z (1,2) = ^ sin20 VZA(l) + (- -> sin20 + — - ~ cos20)VZA(2) p 3 6 6 6 + (- i - sin20 -  +  + (- f- sin20 + —  cos20)VZB(2) + (~„sin20 + ^|)VZB(3)  + (- -|- sin20 +  cos20 + ~)VZC(1) + ( j sin20 - -^|)VZC(2)  1  cos20)VZA(3) + (- |- sin20 -  cos20 - ~)VZB(1)  VJ  + (- - sin20 - -|- cos20)VZC(3)  (1-30)  /2 Z ( l , 3 ) = - j [cos0VZA(l) + cos(0 - 120) VZA(2) + cos(0 + 120) VZA(3) p  + cos0VZB(l) + cos(0 - 120)VZB(2) + cos(0 + 120)VZB(3) + cos0VZC(l) + cos (0 -120)VZC(2) + cos(0 + 120)VZC(3)] Z (2,1) = i sin2@VZA(l) + (- \ sin20 - ^ cos20 - ~ )VZA(2) p J b o b  (1-31)  + (- |- sin2© + ^ ~ cos2 0 + ^)VZA(3) + (- |- sin20 + ^  - ^ ~ cos20)VZB(l) + (- -|.sin20 + ^  + (~ sin20 - -5)VZB(3) + (- \ sin20 o o b + (- | sin20 +  )VZG/2) + (- | sin20 -  cos20)VZB(2)  + ^ cos2o)VZC(l) o b cos20)VZC(3) (1-32)  56.  Z (2,2) = (4 - T cos20)VZA(l) + (-7 P J J + (-7+7 0  c o s 2  0  © +  + (4 + T cos20 + ~ 0  0  + (-7 +  6 6  1 1 ' + 7 cos20 +  0  0  ^ i +  C O s 2 e  ~ "T  + 7 cos20 - ^ D  sin20)VZA(2) O  O  sin20)VZA(3) + (- 7 + T> cos20 6 6  6  sin20)VZB(l)  sin20)VZB(2) + (- 7 - \ cos20)VZB(3) 6 J i i sin20)VZC(l) + ( - 7 - 7 cos20)VZC(2)  O  O  s i n 2 0  )  v z c  j  ( )  (1-33)  3  /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) + s i n (0 + 120)VZB(3) + sinOVZC(l) + sin(0 - 120)VZG(2) + sin(0 + 120)VZC(3)]  (1-34)  /2 Z ( 3 , l ) = - j [cos0(VZA(l) + VZA(2) + VZA(3)) + cos(0 - 120)(VZB(1) + VZB(2) p  + 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)) + s i n (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)] Define E  p  (1-37)  by the following expression: E  p  = [P] V  No  (1-38)  El(l)  [P]v  No  - [P] El(2) El(3) (1-39)  57.  cosG E l ( l ) + cos(0 - 120)E1(2) + cos(0 + 120)E1(3) 3  s i n e . E l ( l ) + sln(0 - 120)El(2) + sin(0 + 120)E1(3)  1_  El(l) +  El(2)  vT  /2  +  El(3) /2 (1-40)  The d-q-0 voltages may now be determined. V  = Z (1,1) i Q  p  d  + Z p  (1,2)1 .  q  + Z (1,3)1  + E (1)  p  o  (1-41)  p  or v  d  =  0j + § cos20)VZA(l) + (- i - i - cos20 + 3  +  - i - -|- cos20 - ~  +  - T - T cos20 + O  + + +  sin20)VZA(3) sin20)VZB(l) D  ~ - j cos20 - ~  sin20)VZB(2) + (-\ + \ cos20)VZB(3)  ~ \ ~ \ cos20 - ~ j - |- cos20 +  + +  D  sin20)VZC(l) + (- |- + -j cos20)VZC(2) sin20)VZC(3)}- i  sin20VZA(l) + (- |- sin20 + - I - sin20 - Sr + 6 6 1  +  sin20)VZA(2)  6  -  d  cos20)VZA(2)  cos20)VZA(3)  /3 /J sin20 - -~ cos2G - — )VZB(1) 6 6  +  - \ sin20 + ^ | cos20)VZB(2) + (^ sin20 + b o 5  +•  _ 1 6  +  -  s i n 2 9  + .if 6  c o s 2 0  +3.  o  )VZB(3)  )VZC(1) + (| sin20 6  3  )vzc(2)  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)}•! 0 [cos0E 1(1) + cos(0 - 120)E1(2) + cos(0 + 120)E1(3)J (1-42) Equation (1-42) can be rewritten as:  • v. = d o  cos20[2VZA(l) - VZA(2) - VZA(3) - VZB(l) - VZB(2) + 2VZB(3)  - 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 f o r V and V gives: q o v  of,  q  •y ^  •h  +  y  2  >  2  )  •\  +  z  p  ( 2  >  3 )  •  h  +  E  ( 2 ) P  (1-44)  59.  =  {~ sin2Gl2VZA(l) - VZA(2) - VZA(3) - V Z B ( l ) - VZB(2)  +  2VZB(3) - V Z C ( l ) + 2VZC(2) - V Z C ( 3 ) ]  /3 +  -g- cos20[-VZA(2) + VZA(3) - V Z B ( l ) + VZB(2) + V Z C ( l ) - V Z C ( 3 ) ]  /3 -  - j [VZA(2) - VZA(3) - V Z B ( l ) + VZB(3) + V Z C ( l ) - V Z C ( 2 ) ] } i  +  {- j COS20I2VZA(1) - VZA(2) - VZA(3) - V Z B ( l ) - VZB(2) + 2VZB(3)  - V Z C ( l ) + 2VZC(2) - VZC(3)]  +  |- [2VZA(1) - VZA(2) - VZA(3) - V Z B ( l ) + 2VZB(2) - VZB(3)  - V Z C ( l ) - VZC(2) + 2VZC(3)]} i q  /2 + ~  {siri0[VAZ(l) + V Z B ( l ) + V Z C ( l ) ]  + sin(0 - 120)[VZA(2) + VZB(2) + VZC(2) + sin(0 + 120)[VZA(3) + VZB(3) + V Z C ( 3 ) ] } i  v  o  = V'* 3  1  '  +  Z  P  (  3  '  2  )  * \  +  Z  P  (  3  '  3  )  * "0  +  E  /2 V Q  =  {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  £  p  ( 3 )  (1-46)  60.  + | (VZA(l) + VZA(2) + VZA(3) + VZB(l) + VZB(2) + VZB(3) + VZC(l) + VZC(2) + VZC(3)} i  + - [ E l ( l ) + El(2) + El(3)] ^  (1-47)  The following alphanumeric names are used i n 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  BCONS2 = —  Z  (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 + s i n ( 0 - 120) • C0NSA2 + sin(0 + 120) • CONSA3 CONSTB = CONIOB • i + ./- (sin0 • E l ( l ) + sin(0 - 120)E1(2) o V J t  + s i n ( 0 + 120)E1(3)) DCONS1 =  (cos0 • CONSA4 + cos(0 - 120) • CONSA5  + cos(0 + 120) • CONSA6) DC0NS2 = —  (sin0 • CONSA4 + sin(9 - 120) • CONSA5  + s i n ( 0 + 120) • CONSA6) Taking advantage of the symmetry i n (1-43) and (1-45), the expressions f o r V,, V and V reduce to: d q o  V  d  = (AC0NS1 + BC0NS1) • i  d  + (AC0NS2 + BCONS2) • i  + CONSTA (1-48)  V  q  = (AC0NS2 - BCONS2) • i , + (-AC0NS1 + BC0NS1) • i + CONSTB d q (1-49)  V = DC0NS1 • i , + DCONS2 • i + CONIAC • i 0 q o d  + - ( E l ( l ) + El(2) + El(3)) /3  (1-50)  J  The p a r t i a l derivatives of » q » ^ v  V  d  0  with respect to the e l e c t r i c angle 9 may  be found quite readily from (1-41), (1-42) and (1-43).  3V • ~  = -2 • AC0NS2 • i , + 2 • AC0NS1 • I  30  d  3V ^-=2 ofc)  - CONSTB  q  • AC0NS1 • i , + 2 • AC0NS2 • i + CONSTA d q  3V -rrp- = - DCONS2 • i , + DC0NS1 • i 30 d q  v.l-31;  (1-52)  n  (1-53)  62.  (1-48), (1-49),  From to  i  d  , i  q  , i  Q  and  (1-50),  the p a r t i a l  derivatives  with  respect  are:  T - r  - =  2  0"54)  AC0NS1 + BC0NS1  d  = AC0NS2 + BC0NS2  -zA o I  "  (1-55)  q  -^A = C0NI0A 9 ,  (1-56)  o  8v ^rA 9  = AC0NS2 -  BC0NS2  (1-57)  d  l  9v = - AC0NS1 + BC0NS1  (1-58)  = CONIOB  (1-59)  = DC0NS1  (1-60)  = DC0NS2  (1-61)  0 = C0NI0C  (1-62)  ^A  oI  q  av •TA a  'o  — ° 8  ,  d  av - r^ ^  ai  q  3v 9 ,  o  APPENDIX  STEADY  I I  STATE PROGRAM  "MSTEA"  64. I TERMINAL  SYSTEM  FORTRAN  G(41336)  MSTEA  SUBROUTINE M S T E A ( I , OMEGA, D E L T A T ) COMMON / V L A O / CM(7,7),CT(3,3),CIF(10),ZTR(70)tZTX{70),ZR<?.0) CCMMON / V L A D / V Z A ( 3 0 ) , V Z R ( 3 0 ) , V Z C ( 3 0 ) , V Z D ( 3 0 ) , E H 3 ) , C U R R T { 3 ) COMMON / V L A O / T H E T < 1 0 J ,ZX<20 I,ZML(90),ZMR(5C),TEM(10),ANG,Z3 COMMON / V L AD/ M H 1 0 ) , M A C H COMMON / M A C H I N / CCNST1,CONST?,CCNST3,HCT,HDTI,ACON COMMON / M A C H I N / A C 0 N 1 , A C 0 N 2 , A C 0 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 I N / R I N E P T 1 8 1 , S T I F F ( 7 ) , D A M p ( 8 ) , TORQUE ( 7 ) , T H E O L D ( 8 ) COMMON / M A C H I N / WOLD I8),KBAND(7,7),HIST1(8),HIST2(8 I COMMON / M A C H I N / NPOLES,NMASS,NROTOR DOUBLE P R E C I S I O N 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 , 1THET,ZX,Z^LfZMR,TEM,ANG.Z3 DOUBLE P R E C I S I O N CONST 1, C 0 N S T 2 , CON ST 3, HOT , H-CT I , ACO N , A C O N 1 . A C 0 N 2 , 3ACON3,AC0N14,ZTY,LPARK,RINERT.STIFF,DAMP,TORQUE,THEOLD,WOLD, 3KBAND,HIST1,HIST2 REAL *8 AS,BS,CS,0S,FS,AR,BR,CF,AI,BI,CI,A,B,C,D,E,F,A1,B1, 3 A L , B L , C l , D L , D E L T A T , OMEGA,T E, A C 0 K 8 , A C 0 N 9 , A C C N 1 0 . A C O N I 1 . A C O N 1 2 , S A C 0 N 1 3 , A C 0 N 1 A , A C O N I B , A CONS 1 , A C O N 2 A , A C O N 2 R , A C C N S 2 , B C D N S 1 , P C O N S T , 3THMECH,WMECH,G(8I,DT(7,7 I ,RZ(8 I.DET.PPCLEStTIME.TORQM INTEGER I P E R M ( 1 4 ) , P A R A M l , P A R A M 2 , P A R A M 3 ISAVE^I N = U*2> /3 WRITE(6,600) N 600 FORMAT ( IH , 2 0 X , • SOLUT ION FOR MACHINE NUMBER ' , 1 4 ) MI(N)=7 IF(ZML(N*9-2) .EQ.O.DO) MHN)=6 MU1 = MI { N ) PARAMl=NMASS-l PARAM2 = N R 0 T 0 R - 1 PCONST=2.DO/DFLOAT(NPOLES) PP0LES=NP0LES/2 IF(N.GT.l) GO TO 6 1 0 HDT=0ELTAT/2.D0 HOT I = 2 . D O / D E L T A T ANG=(?.DO/3.DO)*3.1415 926535^00 CGNST1=DSCRT(3.0C» CCNST2=CCNST1/DSQRT< 2 . 0 0 ) 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  C  ZMR{ I L M ) = ZMR{ I L M | 4 Z R ( I L N ) + 2.D0*ZR( I L N + 1 ) ILM=9*N AS=ZML<ILM)*OMEGA CS=ZVL(ILM-8)*CMEGA DS=ZML(ILM-7) FS=ZML(ILM-4)*0MEGA !LM=5*(N-l) BS = ZMR( ILM + 1 ) CALCULATE P O S I T I V E SEQUENCE CURRENTS ************************ J =I+1 K=I + 2  AR=VZA(I)  BR=VZA(J)  65. TERMINAL  SYSTEM  FORTRAN  MSTEA  G(41336)  CR=VZA(K)  AI=VZB(II  611  BI=VZB(J) CI=VZR(K) AM3R+CR)*0.5D0 B=(C I - B ! > * C O N S T 3 C=(C I4-B I ) * 0 . 5 D 0 D=(BR-CR)*CONST3 BR=(AR*BR*CR)/3.DO BI={AI+RI+CI )/3.DO CR=( A R - , A - e J / 3 . D 0 CI=(AT-C-n)/3.DO AR={AR-A+B)/3.D0 AI=(AI-C*D1/3.00 W R I T E ( 6 , 6 1 l » AR » A I FORMAT{• « , " P O S I T I V E S ' I M A G I N A R Y PART= WRITE(6,612)  612 613  614 :  C  FORM A T ( 1 H  SEQUENCE CURRENT:  REAL  PART=  «,D15.7,5X,  ',015.7)  •'INFORMATION  ABOUT T H E U N B A L A N C E OF T H E S Y S T E M  «)  WRITE(6,613) 8R,BI FORM A T ( ' » , » Z E R O SEQUENCE C U R R E N T : R E A L PA RT = ' . D 1 5 . 7 . 5 X , a»IMAGINARY PART= « , D 1 5 . 7 > WRITE(6,61A) CR,CI FORMAT{• « , ' N E G A T I V E S E Q U E N C E C U R R E N T : R E A L PART=' » , D 1 5 . 7 , 5 X , 3 ' IMAGINARY PART = %D15.7> C A L C U L A T E P O S I T I V E S E Q U E N C E CURRENTS IB A N C IC **************** BR=-0.5D0*AR+C0NST3*AI BI=-C.5D0*AI-CCNST3*AR CR=-0.5 00*AR-CCNST3*4I CI=-0.5DC*AI*CCNST3*AR FIND THE E L E C T R I C A L P O S I T I O N OF THE ROTOR ********************  MN=N*2 CL=ZRIMN-1)-ZR(MN» 0L=ZX(MN-1)-ZX(MN) A L = Z R ( M \ > * ( A R+BR + CR) B L = Z X ( M N ) * ( AI + B I + C I ) AL--AL+BL VZC(I)= V Z C { I » - A L - C L * A R * D L * A I VZC<J)=VZC(J)+AL-CL*BR*DL*BI VZCC K ) = V Z C ( K ) + A L - C L * C R + D L * C I AL=ZR(MNI*{AI*BI*CII BL=ZX(MN)*(AR+BR+CR) AL=AL+BL = vzo( D - A L - C L * A I -OL*AR VZD( J ) = V Z D ( 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=VZDII) C=A*BS*AR-FS*AI 0=B+BS*A I * F S * A R . C=DATAN2(D,Cl E L E C T R I C ANGLE OF T H E R O T O R * * * * * * * * * * T H E T ( N ) = C + 1 . 5707<5632679500 F I N D I D , I C » V D , V Q T C C A L C U L A T E T H E CURRENT HM=N*7-6 A1=DSQRT(A**2+B**2) B i = D A T A N 2 ( B » A) B=B1-C  vzn{ n  C C  IF  *******************  66. TERMINAL  SYSTEM  620  C  C  MSTEA  G(41336>  F=CONST2*Al*DSIN(B) ZTXlILM)=F B=CONST2*Al*OCOS(B) ZTX{ILM+1)=B A1=DSQRT(AR**2+A1**2) Bl=OATAN2(A!«AR) A=B1-C D = C O N S T 2 * A l * D S I N C A) ZTR(ILM)=D A=CCNST2*Al*OCOS(A) ZTR(ILM+1)=A ILM=ILM+2 ZTR(ILM)=0.00 Z T X ( I L M )=O.DO K=MI(N)-3 DO 6 2 0 L = 1 , K J=L+ILM ZTRUI=0.00 ZTX(J)=0.D0 C=(B+8S*A-CS*D)/(DS*0MEGA? ZTX(ILM^1I=-ZMR((N-l>*5*2)*C K=ILM+1 ZTR(K)=C  FORM  62? 621  FORTRAN  INDUCTANCE  MATRIX**********  CALL I N 0 U ( Z M L » C 0 N S T 2 , L P A R K ,N,MU1) ILM=(N-1)*7 DO 6 2 1 J = 1,M|J1 ILN=ILM+J ZTY(TLNI=0.00 DO 6 2 2 K=1,MU1 Z T Y ( I L N ) = Z T Y ( ILN )+LPARK ( J , K J * Z T R (K > CONTINUE  TRANSFORM  VOLTAGES  AND CURRENTS  BACK  TO C O I L  VARI A B L E S * * * * * * * * * *  TE=TMET(N I AC0N8=DC0SfTE) AC0N9=DCCS(TE-ANG) AC0N10=DC0S(TE-ANG) AC0N11=0SIN(TF) ACCN12=DSIN(TE-ANG) AC0N13=DSIN(TE+ANG) AC0NIA=AC0N1*(AC0N8*D+AC0N11*A) A C 0 N I B = A C C N 1 * ( A C C fv 9 * 0 * A C O N 1 2 * A I AC0NS1=ACGN1*(ACON10*D+ACON13*A) AC0N2A=ACCN1*{ACCN8*F+AC0N11*BJ AC0N2B=AC0Nl*(ACCN9*F*ACON12*B) A C 0 N S 2 = A C C N l * ( A C 0 N 1 0 * F +AC0N13*B I • TE = P P 0 L E S * ( Z T Y U L M + 1 > * A - Z T Y ( I L M * 2 » * D ) WRITE<6,623» 0,A,C,F,B,THET<NJ 623 FORM AT I • « , ' I - D = « ,D1 5 .7 , 8 X, • I-C= • , 01 5 . 7 , E X, • I- F= «,D15.7,/, S « V - D = » , 0 1 5 . 7 , 8 X , « V - Q = • , D 1 5 . 7 , 8 X , « E L E C A N G L E (RAD I A N S ) = ' , 0 1 5 . 7 ) C I N I T I A L I Z E THE MECHANICAL E Q U A T I O N S ********** THMECH=PCCNST*THET{N) VMECH=PCCNST*OMEGA WRITE(6,624J TE '624 FORMAT(* • , « E L E C T R C M E C H A N I C A L TCRQUE: ROTOR=•,Dl5.7,5X) IF(NMASS.GT.l) GO TO 6 5 0 B=0.D0 RTNERTt 1 ) = U D 1 0  67. 4 TERMINAL  SYSTEM  FORTRAN  MSTEA  G141336)  OAMP{11*0.00  TnPOUF(l)=TE THEOLO(1)=THMECH WOLDl1)=WMECH G( 1 ) =THMECH GO TO 6 8 2 650  625 630 632 635 626  651  652 655 656  657 658  659 660 661 662  664  ILM=(N-1)*7 TORO M=T E IF I N R O T O R . E Q . N M A S S ) TO ROUE(NROTOR 1 = 0 . 0 0 I F f T O & O U E t N P O T O R ) . E Q . O . D O ) GO TO 6 3 0 TORQUE INROTOR>=-( P P O L E S * C ) * { Z T X l I L M + 4 ) + T O R Q U E T0RQM=TE-T0R0UEINPCT0R) WRITE16,625) TORQUE(NROTOR) F O R M A T l • + « , « E X C I T E R = ' , D 1 5 . 7) DO 6 3 2 1 = 1 , N M 4 S S TORQM=TORQM+WMECH*CAMP(I} DO 6 3 5 I = 1 , P A R A M 2 TORQUE(I)=(TORQUE(11/100.D0)*TORQM WRIT E( 6 , 6 2 6 ) ( 1, TORQUE ( I ) , I = 1 , P A R A M 2 ) FORMAT!'  • , • MECHANICAL  TORQUES  FROM T U R B I N E S * » / » 7 I * ( • , 1 1 » * ) ' »  3D13.6,2X)) IF(NROTOR.EQ.NMASS) GO TO 6 5 2 A=DA MP(NROTOR) B=RINERT(NROTOR) ILN=NR0TCR*1 DO 6 5 1 I = I L N , N M A S S J=I-1 O A M P l J ) = CAMP( I | RINE T(J)=RINFPT{I) CAMPt NMASS) = A RINERTINMASS)=B DO 6 5 5 I = 1 , P A R A M 1 DO 6 5 5 J=l,PARAMl KBANDlI,J)=O.DO DO 6 5 6 I = 1 , P A R A M 1 KBANDII , H = S T I F F ( I > I F f N R 0 T 0 R . L E . 2 1 GO TO 6 5 8 CO 6 5 7 I = 2 , P A R A M 2 J = I- 1 K B A N D l I , I ) = K B A N D ( I , I ) + ST I F F { J ) KBAND(J,I)=-STIFFIJ) K B A N D l I , J ) = K B A N O t J , I> I F I N R 0 T 0 R + 1 . G E . N M A S S ) GO TO 660 I F I N M A S S . E Q . 2 ) GC TO 6 6 0 ILN=NMASS-2 DO 6 5 9 I = NROTOR, ILN J=I*1 KBANDlI,I)=KBANDl1,1)+STIFF(J) KBANDl J , I ) = - S T I F F U ) K B A N D d , J ) = KEAND I J , I ) DO 6 6 1 I = 1 , P A R A M 1 G(II=TOPQUElI)-DAMP(I)*WMECH DO 6 6 2 1=1,NMASS WOLD(I ) = WMECH I F I N R O T O R . E Q . N M A S S I GO TO 6 6 4 GlNROTOR)=GINROTOR)+STIFF INROTOR)*THMECH IFINROTOR.EQ.1) GC TO 6 6 5 I=NR0T0R-1 p  <NROTOR)*C)/WMECH  68. M TERMINAL SYSTEM FORTRAN G H 1 3 3 6 )  MSTEA  6 ( I ) a R ( I J * S T l F F ( I)*THWECH 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 ' , 5 X , 7 ( D 1 5 . 7 ) ) THEOLD(NMASS)=THMECH I F t M R O T O R . F O . 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 665  E = STIFF<n  675  676 677 682  678  679 68C  681  673 674  F=G( I + l ) HIST21I ) = - B * C + ( B + E ) * D - E * F B=E C=D D=F HIST21NMASS)=-B*C*B*D B=HIST2{NROTOR ) I F t N R O T O R . E C N M A S S ) GO TO 677 DO 676 J=NPOTOP,PARAM1 H I S T 2 U J=HIST2I J + l ) UM=(N-1)*7  HIST2lNMASS)=B+PP0LES*tZTY(ILM+l)*ZTR{ILM+2)-ZTY(ILM+2)*ZTR(ILM+1  3) 00 678 J=1,NP0TCP HIST2U)=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) HIST1{I)=(HDTI+DAMP(I)/RINERT(I))*THEOLD{I)+2.D0*WMECH 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) WRITE18) AC0N1 A » A C C N 1 B , A C 0 N S 1 » ACON 2 A, A C 0 N 2 B » A C 0 N S 2 , 3ZTR(ILM+4),TE,A,B,C,THET(N) WRITE(6,673) NVASS,N FORMAT!• »,» I N I T I A L MECHANICAL ANGLES OF MASSES 1 TO ' , 1 1 , • FOR " 3CHINE ' , 1 2 ) WRITE(6,674) { I,G( I),I=1,NMASS) FORM A T ( ' ' , 4 ( • M A S S - • , 1 1 , D I 5 . 7 , 5 X ) 1 I=ISAVE RETURN  APPENDIX  TIME STEP  I II  PROGRAM  1  M TERMINAL  SYSTEM  FORTRAN  SUBROUTINE  0(41336}  GEN  70.  G E N U ,TI M E , D E L T A T , C M E G A ) CM ( 7 , 7 ), CT{ 3 , 3 ),C I F ( 1 0 ) , Z TR ( 7 C ) , Z T X ( 7 0 ) , ZR ( 2 0 ) 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) THET(10),ZX(20),ZML(90),ZMR<50),TEM<10),ANG,Z3 MI(10),MACH  COMMON COMMON COMMON CCMMCN  /VLAO/ /VLAO/ /VLAO/ / V L AD/  COMMON COMMON COMMON COMMON COMMON  / N"ACH IN / /MACHIN/ / M A C H IN/ /MACHIN/ /MACHIN/  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 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 ) 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 ) W O L D ( 8 ) , K B A N D ( 7 , 7 ) , H I S T 1 I 8 ) , H I ST2 ( 8) NPOLES,NMASS,NROTOR  DOUBLE 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 , 1THET,7X,ZML,ZMR,TEM,ANG,Z3 D O U B L E P R E C I S I C N C O N S T I , C O N S 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 , 3AC0N3,ACCN14,ZTY,LPARK,RINERT,STIFF,DAMP,TCRQUE,THEOLD,WOLD, aKBAND,HISTl,HIST2 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 , 3X0LPI7),V0LD(7),PSI0L0(7),X(7),V(7),PSI17),C(15),CCNSA1,CONSA2, 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 , SCONI OA,CONI 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 , A C O N 1 i , A C C N 1 2 , A C 0 N 1 3 , A C O N I A , 3AC0N1B,ACCNS1,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 ) INTEGER IPERM{30),PARAM1,PARAM2,PARAM3 N=(I+2)/3 ISAVE=T MU1=MI(N) T  PARAMUKMASS-1 PARAM?=NP0T0R-1 PARAM3=NMASS+MU1  593  594  595  PP0LES=NPCLES/2 ITER=0 ICONV=0 DO 5 9 3 1 = 1 , N M A S S THENEW{ I)=THEOLD(I)+WOLD(I)*DELTAT !LM={N-1)*7 00 594 1=1,ML1 X O L D U ) = ZTR( ILM+I ) X( I ) = X O L DI I ) VOLD(I)=ZTX fILM+I } PSIOLD(I)=ZTY(ILM+I) DC 5 9 5 1=1,5 ILM=N*5-I+1 ILN=7-I+1 RVEC U L N ) = Z M R { I L V ) ILN=2*N RVEC(2)=RVEC(3)-3.D0*ZR<ILN) RVEC(i)=RVEC(2) I1MSAVE I2=ISAVE+1 13 = 1 S A V E + 2 ,C0NSA1=-AC0N*{VZA(I1)+VZB( I1)+VZC( 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) ) C C N S A 4 = - ( V Z A ( I 1 ) + V Z A ( I 2 ) +V Z A ( 1 3 ) ) C O N S 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) + V Z C I I 2 ) + V Z C U 3 ) ) CONI0C=(1.D0/3 .00 )*(C0NSA4+C0NSA5+C0NSA6) AC0N1A=-2.D0*VZAU1)+VZA( I 2)+VZA(I3)+VZB(11)+VZB(12)-2.D0*VZB(  13  ^ TERMINAL  SYSTEM  FORTRAN  G141336)  71.  GEN  3 + V Z C ( T 1 ) - 2 . D 0 * V Z C ( 1 2 >*VZC( 1 3 ) 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 ) IFfNMASS.EO.lI GO T O 6 0 1 •C***BEGIN C O N S T R U C T I O N O F 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 FIRST*** C * * * U P P E R L E F T SOUARF ( T I M E INDEPENDENT)*** DO 5 9 6 J=1,PARAM1 DO 5 9 7 K = 1 P AR AM I 1  597  JACOB(J,K)=KBANC(J,K)*HOT/RINERT(J)  596 JACOB(J,J)=JAC03(J,J)+HDTI+DAMP(J)/RINERT!J) C * * * U P P E R RIGHT PART (TIME INDEPENDENT)*** DO 5 9 3 J=l,PARANi DO 5 9 8 K=NMASS,PARAM3 598  JACOB(J C=O.DO  t  K) = 0 . D O  B=STIFP(PARAM2) JACOB(PARAM2,PARAM3)=-HOT*B/RINERT(PARAM2) IF(NROTOR.EO.NMASS) GO T O 5 9 9 C=ST I F F ( N R O T O R ) JACOB( NRCT0R,PARAM3 )=-HDT*C/RINERT (NROTOR) C**&LOWER L E F T P A R T , S T I L L TIME INDEPENDENT - EXPLOITS 599 DO 6 0 2 I = N M A S S , P A R AM 3  00 602  602  SYMMETRY***  J=1,PARAU1  JACO^U , J)=JACOB(  J,I)  C***DIAGONAL E L E M E N T ON L O W E R R I G H T 601 JACOBlPARAM3,PARAM3)=HOTI IF(M M A S S . G T . 1 ) J A C O B ( P ARAM 3RINERT{NMASS) C * * * C A L C U L A T I O N OF E L E C T R I C PART 603 THE=PPOLES*THENEW(NMASS)  CORNER  OF  JACOB I A N * * *  3,PARAM3)=HDTI+{CAMP(N^ASS)+HDT*(B+C))/ OF  JACOBIAN  -  LOWER  RIGHT  SQUARE***  WNEWDT=THE-(HOT*WCLC(NMASS)*THEOLD(NMASS))*PPOLES THE2=2.D0*THE AC0N4=ACCN2*0C0S(THE2) , AC0N5=ACCN2*0SIN(THE2) AC0M6=ACCN3*DC0S(THE2) AC0N7=ACCN3*0SIN(THE2) AC0N8=DCGS(THE) ACCM9=DCCS(THF-ANG) ACON10=DCCS(THE+ANG) AC0N11=DSIN(THE) ACCN1?=DSIN(THE-ANG) AC0N13=DSIN(THE+ANG) ACCNS1=ACCN4*ACGMA4-AC0N7*AC0NIB AC0NS2 =A C C N 5 * A C C M A - A C 0 N 6 * A C C N l B BC0NS1=ACCN2*AC0M A-0.5D0*(VZB( I2J-VZRCI 3)-VZC(I 2)+VZC(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 >*VZ3(I 3 >+VZCfI 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 C O N I OB = A C C N ' 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 DCONSl=ACQN*(ACON8*C0NSA4 + ACCN9*CONSA5 +ACON10*CONSA6 ) 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 * C O N S 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)=(ACCNS1*BC0NS1)*X 1 1 ) * ( A C O N S 2 * 8 C 0 N S 2 ) * X I 2 >+CONSTA CONSTB=CONT0B*X(3)•AC0N1*(AC0N11*E1(1)+ACON12*E1(2)+AC0N13*El(3) V(2)=(ACCNS2-BCONS2)*X(1) + (-ACONS1+RCONS11*X(2) + C0NSTR 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 34EK3))  605  V(4)=V0LD(4) DO 6 C 5 J=5,MU1 V(J)=O.DO  N  TERMINAL  SYSTEM  FORTRAN  C***EVALUATE OO 6 0 7  G(41336)  PRESENT  FLUX  72.  GEN  LINKAGES***  J=l,*Ul  PSI(J)=O.DO DO 6 0 6 K = 1 , M U 1 P S I ! J ) = PSI(J)->-LPARK{J,K>*X(KI CONTINUE  6C6 607  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 IF!ICONV.EG.1) GC TO 6 1 5 K=PARAM3-1 C***BEGIM FILLING ENTRIES OF ELECTRIC DO 6 C 8 I = N M A S S , K K1=I-PARAM1 DO 6 0 8 J = N M A S S , K K2=J-PARAM1 6C8 JACOB! I , J )=L P A R K ! K 1 , K 2 J J1 = NMASS J2=NMASS+1 J3=NMASS + 2 K1=PARAM3-1 DO 6 1 0 J = N M A S S , K 1 K2=J-PARAMl  PART  I T I ME  DEPENDENT)***  JACOB(J1,J)=JACOB(Jl,J)*WNEWDT*LPARKI2,K2) 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 ) JACOB(J » J)=JAC OR(J,J)+HOT*PVEC(K2\  610  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 ) JACOB! J l , J 2 ) =JACOB ( J l , J2J +HDT*! AC0NS2+BC0NS2) J A C O B t J l , J3 ) = H D T * C C N I O A JACOB!Jl,PARAM3)=  !(-2.00*AC0NS2*X{1)+2.D0*ACCNS1*X(2)-CONSTB)* HOT +P S I ! 2 ) ) * P P O L E S JACOB 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 ) JACOB!J2,J2) =JACOB(J2,J2)+HDT*! -ACCNS1+BCCNSI) 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  3  JACOB!J3, Jl>=H0T*CCCNS1 JACOB!J3,J2)=H0T*DC0NS2 J A C O B ! J 3 , J 3 ) = J A C C B ( J 3 , J 3 ) + HD T * C ON I OC  611  J A C O B ! J 3 ,PARAM3) = ( - 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=NMASS+3 DO 6 1 1 I = J , K 1 JACOB!!,PARAM3)=C.D0  C * * * C A L C U L A T E T I M E D E P E N D E N T ROW U N D E R E L E C T R I C BLOCK IF(NMASS.GT.l) GC T O 6 0 9 RINERTl 1)=1.010 D A M P ( 1 ) = 0 .DO 609 0 0 6 1 2 1=1,MU1 612 JAC0B!PARAM3,PARAM1+I)=LPARK!1, I ) * X ! 2 ) - L P A R K ! 2 , 1 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 )  OF  JACOBIAN***  )*X!I)  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 , M U 1 61A  JACOB I P A R A " 3 ,PAR AMI*I )=JACOB (PARAM3, P A R A M l ^ I  >*HDT*PPOLES/RINERT  3!NMA S S ) . C B E G I N C A L C U L A T I O N O F R I G H T HAND C***MECHANICAL PART FIRST*** 615 IFINROTOR.EO.l) GO T O 6 1 7 DO 6 1 6 I = 1 , P A R A M 2  SIDE  VECTOR  •G«  * * * * *  N TERMINAL 616 617  618 619  SYSTEM  FORTRAN  73.  G(41336)  GEN  T H E T A S I I ) = T H E N E W ( I) T H E T A S ( N R O T O R ) = T H E NEW(NMASS) IF(N^ASS.EO.1) GO T O 6 2 2 I F ( N R O T O R . E C . N M A S S ) GO T O 6 1 9 K=NR.OTOP-l DO 6 1 8 I = K , N M A S S THETAS( I } =THENEW(1-1) D=THPTAS(1) DO 6 2 0 I = 1 , P A R A M 1 E = ST I F F ( I ) F=THETAS(1 +1 )  620  HNEW2( I ) = - 8 * C + ( B - * E ) * 0 - E * F B=E C=D D=F HNEW21NMA5S)=-B*C+B*D B=HNEW2(NROTOR ) IF(NROTOR.EC.NMASS) GO T O DO 6 2 1 J = N R 0 T 0 R , P A R A M 1  621 622  622  HNEW2(J)=HNEW?(J+l) HNEW2(NMASS)=B*TE DO 6 2 3 J = l , N R O T O R 623 HNEW2 ( J ) = H N E W 2 ( J ) - T O R Q ' J E ( J ) DO 6 2 4 1=1,NMASS H N E W 2 ( I ) = H N E W 2 I I ) * H D T / R I N E R T (I ) B=(HDTI+CAMP{Il/RINERT tI))*TrENEW(I ) HNEW1 ( I ) = B C=B-HIST 1{I) 624 G(I)=-C-HNFW?(I)-HIST2(I) G(PARAM3)=G(NMASS) IF(ICONV.EQ.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=PARAM3-1 DO 6 2 5 J = N M A S S , L K=J-PARAM1 625 G(J)--PSnK)+PSIOLD(K)-HDT*(V(K)+VOLD(K)+RVEC(K)*(X(K)+XOLD(K))) C=THENEW(J1)-THEClC(Jl> 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)-PS I OLD(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 I O L D ( 11 C * * * S O L V E THE S E T OF L I N E A R EQUATICNS (UP T O 15 X 15)*** CET=1.D-14 IJK = 0 CALL 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 T O 6 3 1 DO 6 3 3 1=1,PARAM3 633 IF(DABS(DELTAX(I)).GT.l.D-9) IJK=IJK+i I F ( N M A S S . E Q . l ) GC TO 6 3 7 DO 6 3 8 J=1,PARAM1 638 THENEVM J ) = T H E N E W ( J ) + D E L T A X ( J ) 637 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 ) ,DO 639 J=1,MU1 K1=PARAM1+J 639 X( J ) =X( J ) + D E L T A X ( K l ) ITER = ITER + 1 IF(ITER.GT.10) WRITE(6,656) 656 F O R M A T { • • , » * * * * * S O L U I C N NOT W I T H I N E P S I L O N A F T E R 10 ITERATIONS. 3LAST VALUES USED*****') T  IF(ITER.GT.10)  GO T O 6 4 1  N TERMINAL  SYSTEM  FORTRAN  IF( I J K . G T . O ) WR!TE(6 ,640I  0(41336) GO T O ITER  GEN  603  640  FORM A T ( • ' , • S O L U T I O N BY NEWTCN* ' S METHOD C O N V E R G E D I N ' , I 2 , IT E f (DTICNS • ) 641 IC0NV=1 GO T O 6 0 3 642 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 ) 645 H!ST2(I)=HNEW2(I ) 1 L M = ( N - l )*7 DO 6 4 4 J=1,MU1 ZTR(ILM+J) = X(J ) ZTX( ILM + J ) = V ( J ) 644 ZTY(ILM+J)=PSI(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 T O P H A S E DOMAIN W=AC0\'1*ACCN14*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 AC0N1B=ACCN1*(ACCN9*X(1)+AC0N12*X(2))+W A C 0 N S 1 = A C C N 1 * ( A C C M O * X ( l > ~ A C 0 N 1 3 * X ( 2 ) )*W W=AC0N1*«C0N14*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 AC0NS2=ACCN1*<ACCN10*V(1)+AC0N13*V{2) ) +W T E = P P Q L E S * ( P S I ( 1 > * X ( 2 ) - P S I (2 ) * X ( 1 ) I CURRT(1 )=-ACONlA CURRT(2)=-ACONlB CURRT(3)=-ACCNSl WRITE(6,650) AC0N1A,X(1),X(3),X(4),TE,WOLD(NMASS ) 65 0 F O R M A T ( • • , M - A = • , C I 4 . 7 , 2X , « 1 - 0 = • , DI 4 . 7 , 2 X , • I - 0 = » , DI 4 . 7 , 2 X , « I - F = aD14.7,2X,«T-E= ,014.7,2X,«R0T0R SPEED=«,014.7) i T H E T ( N ) = P P O L E S * T r - ENEW ( NM A S S ) - 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 IF(NMASS.FO.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 ) )* S T I F F < N R Q T O R - 2 ) IF(NMASS.GT.NROTOR) D=(THETAS(NROTORI-THETAS(NROTOR+  aoR) 651  WRITE<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 , aX(4),TE,B,C,D,THET(N) GO T O 6 5 5 631 WRITE(6,653) 653 F O R M A T ! • « , ' S O L U T I C N F A I L E D , CHECK YOUR MATRICES') 655 CONTINUE I=ISAVE RETURN END r i C N S IN E F F E C T * 10,E8CDIC,SOURCE,NOL1ST,NODECK,LOAD,NOMAP ("IONS IN E F F E C T * NAME = G E N , L I N E C NT = 60 4TISTICS* SOURCE STATEMENTS = 2 4 9 , PROGRAM S I Z E = 13032 VTTSTICS* NO D I A G N O S T I C S GENERATED *S IN GEN  1))*STIFF(NRC  

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}]}"
                            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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065483/manifest

Comment

Related Items