THE PROBLEM OF FREQUENCY DEPENDENCE IN TRANSMISSION LINE MODELLING by JOSE R. MARTI E l e c . Engr., Cent r a l U n i v e r s i t y of.Venezuela, 1971 M.E.E.P.E., Rensselaer P o l y t e c h n i c I n s t i t u t e , 1974 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n THE FACULTY OF GRADUATE STUDIES Department of E l e c t r i c a l Engineering We accept t h i s t h e s i s as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA A p r i l 1981 © Jose R. M a r t i , 1981 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by h i s or her representatives. I t i s understood that copying or p u b l i c a t i o n of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 DE-6 (2/79) ABSTRACT In this work, the accurate representation of transmission lines for the d i g i t a l simulation of electromagnetic transients in power systems has been examined. A model has been developed that accounts for the frequency dependence and distributed nature of the line parameters over the entire frequency range. This model can easily be incorporated into a time-domain network solution of the complete power system. The model consists simply of a constant resistence in parallel with a current source evaluated at each time step of the solution. The equivalent resistance results from a f in i t e-step-width discretization of the differential equations of a resistance-capacitance (R-C) network that simulates the line characteristic impedance. The equivalent current source accounts for the time delays and attenuations of the different frequency components of the travelling waves and for the discretization of the time-domain equations. Rational-function approximations are used to synthesize the R-C network and the line propagation ("weighting") function in the frequency domain. These rational approximations allow the corresponding time-domain functions to be obtained directly in a closed-form, thus circumventing the need for numerical inverse Fourier transformations. The numerical technique used to obtain the rational functions yields very accurate, high-order approximations. This technique is based on a direct, step-by-step allocation (and reallocation) of poles and zeros and avoids the instability problems which can be encountered with optimization techniques based on search methods. i i A series of analytical evaluations and simulation tests were performed in order to assess the va l i d i t y of the model. The results of these tests show that the model is accurate, fast, and reliable. The model was incorporated into the code of the University of Brit i s h Columbia's version of Dr. H.W. Dommel's Electromagnetic Transients Program (EMTP). i i i TABLE OF CONTENTS ABSTRACT ACKNOWLEDGEMENT INTRODUCTION General P e r s p e c t i v e on the Problem of Transient S t u d i e s , 1 Purposes and Objectives of the Present Work, 5 CHAPTER 1: Wave Propagation i n Transmission Lines 1.1 General Line Equations, 6 1.2 S o l u t i o n of the Line Equations i n the Frequency Domain, 11 1.3 S o l u t i o n of the Line Equations i n the Time Domain, 13 CHAPTER 2: Simu l a t i o n of Electromagnetic Transients 2.1 Complete Network S o l u t i o n : Time Domain v s . Frequency Domain, 19 2.2 Network S o l u t i o n i n the Time Domain, 20 CHAPTER 3: The Problem of Frequency Dependence 3.1 General Considerations, 24 3.2 Admittances Formulation: Budner, 25 3.3 T r a v e l l i n g Functions Approach: Snelson; Meyer and Dommel, 30 3.4 New Formulation, 37 CHAPTER 4: Synthesis of the C h a r a c t e r i s t i c Impedance and Weighting Function 4.1 General Considerations, 50 4.2 Synthesis of Equivalent Network to Represent the C h a r a c t e r i s t i c Impedance, 52 i v A.2.1 General Considerations, 52 4.2.2 Rational Approximation of Z ^ C C J ) , 56 4.2.3 Equivalent R-C Network, 57 4.3 Weighting Function and History Convolution, 59 4.3.1 General Considerations, 59 4.3.2 Recursive Convolution, 61 4.3.3 Rational Approximation of A^(w), 65 4.4 Numerical Synthesis of Rational Approximations for Z c (C J ) and A ^ ( a ) ) , 66 4.4.1 General Considerations, 66 4.4.2 Approximation Techniques by Other Authors, 69 4.4.3 Rational Approximation Technique in the New Formulation, 70 CHAPTER 5: Implementation of the Solution 75 5.1 Recapitulation, 75 5.2 Evaluation of the Weighted History Source, 76 5.3 Equivalent Impedance, 78 5.4 Some General Comments on Numerical Integration, 82 5.5 Reduced Equivalent Circuit, 83 5.6 Algorithm of Solution, 85 5.7 I n i t i a l Conditions, 86 5.8 Interface with the EMTP, 89 CHAPTER 6: Numerical Results 93 6.1 Introduction, 93 6.2 Frequency-Dependent Line Parameters, 93 6.3 Simulation of the Characteristic Impedance, 95 6.4 Approximation of the Weighting Function, 98 6.5 Analytical Tests: Frequency Domain, 103 6.6 Analytical Tests: Time Domain, 109 6.7 Examples of Transient Simulations, 114 6.8 Field Test Simulation, 119 List of Graphs, 122 CONCLUSIONS 187 APPENDICES APPENDIX I: Some Matrix Relations 189 1.1 Diagonalization of the Propagation Equations, 189 1.2 Diagonalization of the Drop Equations, 190 APPENDIX II: Integration Coefficients 194 11.1 Recursive Convolution Property, 194 11.2 Dommel's Model for R-C Network, 196 BIBLIOGRAPHY 1*58 v i ACKNOWLEDGEMENT For the completion of this work, I would like to say thanks, To Gail, who worked at i t as hard as (or harder than) the author. To Mae, who tried to teach me the idiosyncrasies and beauty of the English language, but who is not to blame for the pupil's limitations. To my brother, Luis, who shares a dream with me. To Dr. Dommel, whom I admire. To Dr. Moore and "his g i r l s " , for making me feel not a stranger in the Department. To the people responsible for the computing f a c i l i t i e s at the University of British Columbia, who make computers a "bit " less mysterious. To Central University of Venezuela who financed my leave of absence at U.B.C., for their support. To my parents. v u INTRODUCTION GENERAL PERSPECTIVE ON THE PROBLEM OF TRANSIENT STUDIES The size and cost of the equipment in a power system is largely de-pendent on the stresses that occur during transient periods. The equipment must have enough insulation to withstand the overvoltages, enough thermic capacity to withstand the overcurrents, and enough mechanical strength to withstand the increased mechanical forces. Transients occur in a power system when the normal steady-state operation is disturbed. The state of the system can be disturbed by nec-essary internal switching operations, such as those needed for connecting or disconnecting the various system components, or by undesirable extrinsic per-turbations, such as atmospheric discharges, faulty insulation, and mechanical breakdowns. Due to the inertia of i t s components, the system cannot, after a perturbation, instantly adapt i t s e l f to the new^conditions. Depending on the nature of the disturbance and the parameters of the system, the transient processes involved can be very fast, with the currents and voltages reaching very high magnitude peaks and having "arbitrary" waveforms. A series of protective devices, such as lightning arresters and fault-detection relays, are normally employed to limit the magnitude and dura-tion of transient conditions produced by undesirable perturbations. These devices can be quite effective in protecting the system, but at the same time, their location and adjustment can be quite c r i t i c a l . A premature or untimely action on their part can cause interference with switching operations and thus jeopardize the normal functioning of the system. The overvoltages produced 1 during switching operations can also, to a certain degree, be controlled. A common practice for this purpose i s the pre-insertion of resistors before connecting or disconnecting a given part of the system. The effectiveness of protective strategies in moderating transient conditions and the ab i l i t y of the equipment to sustain the "remaining" trans-ients can only be properly assessed i f the behaviour of the system during these periods is accurately described and predicted. Tolerances in view of uncertainties can represent considerable increases in equipment costs with no guarantee of optimum operation. This becomes even more important for very high voltage systems where the dimensions of the equipment and their associated financial and environmental costs are very c r i t i c a l . Electromagnetic transients in power systems can be very fast and peak values may occur within microseconds. When the periods of the signals are comparable to or smaller than the t r a v e l l i n g times of the electromagnetic waves, as i t is usually the case in transmission lines, lumped-parameter re-presentations are no longer adequate and the mathematical description of the phenomena must be formulated in terms of partial differential equations in time and space. In transmission lines the conductors are coupled to each other and, therefore, the formulation results in a system of mutually coupled partial differential equations. For purposes of analysis, each component of an n-phase system is normally represented by an equivalent mathematical model of order n. For transmission lines with ground return, the parameters of these models may be very sensitive to variations in the frequency. Since during transient conditions the signals usually contain a wide range of frequencies, the 3 accurate description of the phenomena involves a system of coupled partial differential equations with non-constant coefficients. Due to the complexity of the problem, early studies of transients in power system (before the advent of di g i t a l computers) had to be restricted to very small systems, and involved a large number of simplifications. (Despite these limitations, some of these early works (e.g. [1] to [5]) are of great merit, and many of their practical results are s t i l l being used in industry today.) The need for more general and accurate solutions led to the development of analogical representations in the form of Transient Net-work Analyzers (TNA) (Peterson, 1939, [6]). In these representations the system components are simulated by bench-size physical models. TNA's represented a large step forward from manual computations and made possible the study of more complex problems with less simplifications. As systems grew in complexity and size, however, TNA's became more elaborate and cumbersome. As a result, i t s use became more restricted to specialized centers which could provide the required physical f a c i l i t i e s and trained personnel. (Some of these centers are s t i l l in operation today (e.g. [7], p.3).) In addition to the size and cost factors, TNA's are also limited in their a b i l i t y to simulate real physical systems. One of these limitations is that distributed-parameter elements, such as transmission lines, are re-presented by lumped-parameter models. To diminish the error, the line is re-presented as segmented into smaller sections, but a very large number of sections would be necessary for an accurate representation of the high frequency components. Another main limitation of TNA's is the modelling of frequency dependence. 4 The most important advances i n s o l v i n g the t r a n s i e n t s i m u l a t i o n problem are probably the ones made i n the l a s t two decades. The development of d i g i t a l computers, numerical methods of a n a l y s i s , and e f f i c i e n t p rocessing algorithms have made p o s s i b l e the study of l a r g e r and more complex systems w i t h l e s s s i m p l i f i c a t i o n s . Probably the most important c o n t r i b u t i o n s i n t h i s new period have been the works by Wedepohl (1963, [ 8 ] ) , Hedman (1965, [ 9 ] ) , and Dommel (1969, [10]). Wedepohl and Hedman (independently from each other) a p p l i e d the mathematical t h e o r i e s of matrix manipulation and modal decoupling to the a n a l y s i s of multiphase systems. These methods a l l o w the r e d u c t i o n of coupled multiconductor systems i n t o equivalent single-phase u n i t s . Dommel developed a general computer program capable of s i m u l a t i n g a large v a r i e t y of t r a n s i e n t c o n d i t i o n s f o r systems of any s i z e and complexity. In t h i s program the d i f f e r e n t i a l equations d e s c r i b i n g the t r a n s i e n t behaviour of the system components are transformed i n t o d i s c r e t e models a t f i n i t e time increments. A complete network s o l u t i o n i s then formulated to solve f o r the system c o n d i t i o n s at each time step. The problem of frequency dependence has been the subject of much study i n the l a s t decade. Some of the most important c o n t r i b u t i o n s to the s o l u t i o n of t h i s problem are l i s t e d i n references [11] to [17]. Despite the progress made i n t h i s area, the formulations implemented at the present time s t i l l present a s e r i e s of numerical problems, and t h e i r use re q u i r e s many p a r t i c u l a r c o n s i d e r a t i o n s . These f a c t o r s have l i m i t e d the usefulness of these methods f o r a general c l a s s of t r a n s i e n t s i m u l a t i o n s . 5 PURPOSES AND OBJECTIVES OF THE PRESENT WORK It has been the purpose of the work presented in this Dissertation to try to develop a formulation for the problem of frequency dependence of transmission systems that is accurate, reliable, and easy to apply for a general class of transient simulations. In the development of this work the previous contributions to the solution of the problem have been studied, and their advantages and disadvantages taken into consideration. To achieve the mentioned objectives, the following general philoso-phical precepts have been observed: i) The mathematical complexity of the problem should not obscure the physical concepts involved, i i ) Within the boundaries of the problem studied, the formulation should be of a general nature; that i s , i t should not be restricted to particular cases or conditions, i i i ) The numerical methods of solution and computer routines should have the following characteristics: a) The computer routines should be easy to implement and interface with existing general transient programs (namely, those based on Dommel's formulation). b) The numerical methods should be fast and accurate, but, increased accuracy should be limited to the extent where the additional time, cost, and effort involved are j u s t i -fied by corresponding improvements in the f i n a l results. c) The general nature of the computer routines should be easy to understand, and their external characteristics easy to manipulate by the users of the transient programs. CHAPTER 1 WAVE PROPAGATION IN TRANSMISSION LINES 1.1 General Line Equations A transmission line in an n-phase power system may have a certain number of conductors per phase, ground wires, and ground return (fig- 1.1(a))-However, i f the line equations are formulated in the frequency domain, the 3Sn w a 9 n (a) (b) Fig. 1.1: Transmission line. (a) Original system. (b) Equivalent system. original system can be reduced by matrix manipulation techniques to a system with only one equivalent conductor per phase and a neutral return ( f i g . 1.1(b)). The effects of the bundling of phase conductors, the presence of ground wires, and the behaviour of the ground return path are a l l included in the equivalent impedances and admittances of the reduced system. The line equations, in terms of the reduced equivalent system, are reviewed next. To simplify the notation, capital letters are used to represent 6 7 matrices in the frequency domain; for example, Z a a ( w ) Z a b ( u ) Z a c ( w ) ha™ Z b b ( w ) Z b c ( u ) z c a M z c b ( . ) z c c ( u > ) V ph V a ( w ) V b ( U ) and I ph I a ( U ) I b(.) Ic(«) Voltages are measured from phase to neutral, and impedances and admittances are expressed in per unit length. The relation between voltage and current along the line is given by dV and dx ph ph Z I (series voltage drop equation) , (1.1) dl ph Y V , (shunt current drop equation) dx ph ph Differentiating with respect to x, the wave propagation equations in the frequency domain are obtained: (1.2) and d2V dx 2^ = (Zph'X'ph^ V ^ (voltage propagation equation) , d 2I dx = (Y Z , ) I , (current propagation equation), z ph ph ph (1.3) (1.4) Since a l l the quantities are matrices, these equations constitute systems of coupled differential equations. The equations can be decoupled, however, through matrix transformations(e.g. Wedepohl [8] and Hedman [9]). Some of the 8 properties of these transformations are reviewed next. From matrix theory, a square matrix A can be diagonalized (under certain conditions) by a transformation of the form T _ 1 A T , (1.5) where T is obtained from the eigenvectors of A. Although the s t r i c t mathematical conditions under which the diagonalization is possible can be elaborate, the product ( z ph Y p h) in equation 1.3 is usually diagonalizable. Assuming this is the case, i t is shown in Appendix 1.1 that the inverse product ^) in equation 1.4 is also diagonalizable. That i s , i f P = diagonalizing matrix for (Z , Y , ) , (1.6) pn ph then Q = diagonalizing matrix for (Y ,Z , ) (1-7) exists. A general relationship between P and Q (App. 1.1) is given by Q = ( P V 1 D c , (1.8) where (P*") ^ is the inverse of P transposed, and is an arbitrary, constant, diagonal matrix. To introduce the diagonalized forms into the propagation equations (1.3 and 1.4), a new set of variables (modal quantities) are defined by similarity (modal) transformations: V = P V ( 1 . 9 ) ph m and Tn h = Q zm • ( 1 - 1 0 ) pn m 9 Substitution of these transformations into the propagation equations gives — = [P~ (Z Y )P]V = D V (1.11) ^ ^ ph ph m zy m 7 and d 2 1 [ Q " ( Yph Zph ) Q ] Im= Dyz Xm > ' where, as shown in Appendix 1.1, D = D . (1.13) zy yz With the change of variables indicated in 1.9 and 1.10, the drop equations (1.1 and 1.2) become dVm ,„-l dx V i ph x / m and = (P Z . Q) I (1.14) it - ( Q _ 1 Yph P ) Vm ' ^' 1 5> It is shown in Appendix 1.2 that the coefficient matrices in these equations are diagonal, that is P 1 Z p h Q = Dz (diagonal) (1.16) and Q 1 Y p h P = Dy (diagonal) . (1.17) 10 Matrices and represent the equivalent series impedances and shunt admit-tances of the decoupled modal system: D = Z (modal series impedances) (1.18) z m D = Y (modal shunt admittances) (1.19) y m Matrix 1.13 gives the propagation constant for the different line modes: 2 D = D = Y (propagation constant squared) (1.20) zy yz With the indicated properties of the modal transformations, the line equations in the mode-frequency domain can be written as: dV —r- = Z I (voltage drop equation), (1.21) d x m m dl Y V (current drop equation), (1.22) dx mm d 2v ? m 2 2 - Y V m (voltage propagation equation), (1.23) dx d 2 i m 2 2 = Y 1 (current propagation equation). (1.24) dx A l l these equations are decoupled. Differentiating equations 1.21 and 1.22 and comparing with equations 1.23 and 1.24, the propagation constant Y c a n D e obtained directly as Y = Jz7T~ . (1.25) m m To recapitulate: As described by equations 1.21 to 1.24 the original coupled system of figure 1.1 has been converted by the modal transformations 1.9 and 1.10 into a new decoupled system. Each one of the component modes of this decoupled system can be studied individually as a single-phase system. 11 Accordingly, in what follows in this work, the line solutions are formulated for single-phase lines, with the understanding that these solutions apply to any of the decoupled modes. 1.2 Solution of the Line Equations in the Frequency Domain With a single-phase line as the reference unit ( f i g . 1.2), the solu-tion to the line equations in the frequency domain is next formulated. (Lower case letters are used to denote quantities in the time domain whereas capital letters are used to denote quantities in the frequency domain.) X — i k ( t ) i m ( t ) ^ _ — Ik(«) I m M — k o O rn k O om f t. » t v k (t) v m ( t ) V k M V m N by and (a) (b) Fig. 1.2: Single-phase line (or component mode). (a) Time domain, (b) Frequency domain. In figure 1.2(b) the voltage and current drop equations are given g-Z'I (1.26) 4^ = Y'V, (1.27) dx where and 12 Z' = R' (u) + ju)L'(co) (R', L 1 = series resistance and inductance per unit length), Y 1 = G'(u)) + jwc'(a)) (G', C' = shunt conductance and capacitance per unit length), V = V(u>,x) (line to neutral voltage), I = I(u,x) (series current). (It is emphasized that in the most general case a l l the parameters are fre-quency dependent and the variables are functions of frequency and space.) Differentiating equations 1.26 and 1.27, 2 ^-j = (Z'Y') V (1.28) dx and ^-j = (Y'Z') I . (1.29) dx From these equations, and with 1.26 and 1.27, the general solutions for voltage and current are given by V(oo,x) = A e ^ ( a ) ) x + B e Y ( w ) x (1.30) and T/ \ A -y(w)x B y(w)x . . I ( u ' x ) = z ~ u y e - " M ^ 6 ' ( 1 - 3 1 ) c c where Y(W) = /z ' (w) Y' (co) = a(w) + j3(io) = propagation constant, (1.32) a ( g ) = attenuation, 3(w) = phase displacement , and Z (ID) = = characteristic impedance. (1.33) c / Y (co) The constants A and B in equations 1.30 and 1.31 are defined by the boundary conditions at any point of the line. At the receiving end (x = I), 13 V = A e" y £ + B e y £ m and from where, and - I - A _ e - Y * - '5- > m Z Z c c V - Z I , , m c m. A = (: 2 ) e V + Z I B = ( m 2 c m) As a function of the conditions at node m, the voltage and current at node k are then given by yl -yl yl -yl \ = VJ~—T1—> - Z, T J - — i - 2 — > k m 2. c m 2 and V yl -yl yl , -yl k Z k 2 ; V 2 ; , c ' or in terms of hyperbolic functions (and emphasizing the dependence on tu) , Vk(co) = cosh[Y(co)£]Vm(u)) - Zc(OJ) sinh[y(oi)£] I (oi) (1.34) and = 7~7~T sinh[y(w)£] V (w) - cosh[Y(o)) *] I (w) . (1.35) K. L \m) m m c These equations represent the "cl a s s i c a l " solution to the line equations in the frequency domain (e.g. Woodruff [18], pp. 93-120). 1.3 Solution of the Line Equations in the Time Domain The d i f f i c u l t y i n obtaining a complete solution to the line equations in the time domain can be illustrated from their frequency domain form. Considering, for instance, equation 1.28 14 = (R' + j u L ' X G ' + jo)C')V dx = [(R'G' - (A ' C ) + j u ( R ' C ' + G'L')]V, (1.36) and n o t i n g t h a t i n g e n e r a l the parameters a r e a c e r t a i n f u n c t i o n of f r e q u e n c y , i t can r e a d i l y be seen t h a t o b t a i n i n g the c o r r e s p o n d i n g time-domain form i s not s t r a i g h t f o r w a r d . B e f o r e p r o c e e d i n g w i t h t h i s p o i n t , i t may be u s e f u l to b r i e f l y examine the n a t u r e of the f r e q u e n c y dependence of the parameters i n a r e a l system. In the most elementary case of a c o n d u c t o r c a r r y i n g an e l e c t r i c a l c u r r e n t , the d i s t r i b u t i o n of t h i s c u r r e n t a c r o s s the c o n d u c t o r i s not u n i f o r m , but depends on the r a t e o f change ( f r e q u e n c y ) of the c u r r e n t . The h i g h e r the f r e q u e n c y the more the c u r r e n t i s d i s p l a c e d towards the o u t s i d e of the con-d u c t o r , thus changing the r e s i s t a n c e and i n t e r n a l i n d u c t a n c e ( s k i n e f f e c t ) . The a n a l y t i c a l s o l u t i o n s d e s c r i b i n g the v a r i a t i o n o f these parameters w i t h f r e q u e n c y are f a i r l y c o m p l i c a t e d and i n v o l v e B e s s e l f u n c t i o n s (e.g. Woodruff [18] , pp. 53-710). N e v e r t h e l e s s , i n t h i s elementary case the most i m p o r t a n t component of the t o t a l impedance i s the e x t e r n a l i n d u c t a n c e (which i s not a f f e c t e d by the f r e q u e n c y ) and s t u d i e s t h a t n e g l e c t f r e q u e n c y dependence can s t i l l g i v e a c c e p t a b l e r e s u l t s . A more c r i t i c a l case i n power t r a n s m i s s i o n systems i s t h a t of a phase c o n d u c t o r w i t h ground r e t u r n . In the case of the ground, which can be imagined as a huge co n d u c t o r , the r e s i s t a n c e and i n t e r n a l i n d u c t a n c e con-s t i t u t e a l a r g e p o r t i o n o f the t o t a l impedance, and the e f f e c t o f the f r e -quency i s thus v e r y i m p o r t a n t . The e x p r e s s i o n s f o r the v a r i a t i o n of these parameters w i t h f r e q u e n c y a r e a g a i n r a t h e r c o m p l i c a t e d and i n v o l v e B e s s e l f u n c t i o n s . (The s o l u t i o n f o r t h i s case was f i r s t g i v e n by Carson, 1926 [ 1 9 ] , and, i n d e p e n d e n t l y , by P o l l a c z e k , 1926 [ 2 0 ] . ) . 15 Going back to equation 1.36, the complicated form of the functions R' (to) and L' (co) makes direct transformation of this equation into the time domain very d i f f i c u l t , i f not practically impossible. (The parameters G' and C' can be considered for practical purposes as independent of the frequency.) Despite not being able to obtain a closed-form solution,it i s s t i l l interesting, for the purpose of understanding the phenomenon, to consider the form of the line equations in the time domain under simplifying conditions. To inverse transform the frequency-domain equations, the following property is f i r s t recalled: d nf(t) V n v ( -•< • (jw) F(w) . dt The following conditions are next considered: a) With R', G1, L' and C' independent of to, the inverse transform of equation 1.36 is given by 2 2 = (R'G')v + L'C' + (R'C1 + G'L') . (1.37) 9x 2 9 t 2 3 t A similar expression is obtained from equation 1.29 for the current. From equation 1.37 i t can be seen that even with the simplifying assumption of frequency-independent parameters, the resulting time domain form is s t i l l complicated. b) Neglecting the losses (that i s , R' = G' =0), the inverse trans-form of equations 1.28 and 1.29 is now given by 2 2 ^ - L ' C - ^ - f (1.38) 3x at and - ^ | = C ' L ' - ^ i . (1.39) 3x 3t In a similar way, for the drop equations (1.26 and 1.27), 16 8x " L ( 1 - 4 0 ) and E q u a t i o n s 1.38 to 1.41 a r e the " c l a s s i c a l " e q u a t i o n s f o r wave p r o p a -g a t i o n i n a l o s s l e s s l i n e . They c o u l d have been f o r m u l a t e d d i r e c t l y i n the time domain, but were i n s t e a d d e r i v e d from the f r e q u e n c y domain i n o r d e r to emphasize the l i m i t a t i o n s under which they a p p l y . A g e n e r a l s o l u t i o n to these e q u a t i o n s was f i r s t g i v e n by d'Alembert i n 1747 i n the form and v = f ^ x - at) + f ( x + a t ) (1.42) 1 = BT f l ( x " a t ) ~ R~~ f 2 ( x + a t ) ' (1.43) c c where = / c R / = surge impedance ( r e s i s t a n c e ) and a = J[Tjjr = v e l o c i t y of p r o p a g a t i o n The f u n c t i o n s f ^ and f ^ a r e d e f i n e d from boundary and i n i t i a l c o n d i -t i o n s . Combining e q u a t i o n s 1.42 and 1.43 the f o l l o w i n g r e l a t i o n s h i p i s ob-t a i n e d : v + R £ i = 2 f 1 (x - a t ) . (1.44) From t h i s r e l a t i o n , i t can be seen t h a t f o r (x - a t ) c o n s t a n t t h e q u a n t i t y (v + R c i ) i s a l s o c o n s t a n t . The q u a n t i t y (v + R £ i ) a t p o i n t 1 on the l i n e w i l l "appear" a t p o i n t 2 (at a d i s t a n c e Ax from p o i n t 1) a f t e r a time 17 This can be interpreted by saying that the quantity (v + R^i) "travels" along the line with a velocity "a". The time needed by this "wave" to travel from one extreme of the line to the other is given by Relation 1.44 is the basis for Bergeron's "Method of Characteristics", which is a graphical approach to relate currents and voltages at different points on the line. that the present value of the quantities at node k (fi g . 1.2(a)) is related to the past value of the quantities at node m by That i s , the relation between current and voltage at node k at a given time t is known i f the past values of the corresponding quantities at node m at T units of time earlier have been recorded. With this relation, the equivalent circuits shown in f i g . 1.3 are obtained. (Similar representations are obtained for node m.) To take the losses into account Dommel represents them as T = line travelling time line length a (1.45) From the above interpretation of equation 1.44, Dommel [10] observed v (t) + R (- i (t)) = v (t - x) + R i (t - T) . K. c K m c m (1.46) O (a) (b) Fig. 1.3: Dommel's representation of a transmission line as seen from node k at time t. (a) Series form. (b) Shunt form. (No frequency dependence or distributed losses are considered.) 18 concentrated at the ends and middle of the line. No frequency dependence is considered in this basic representation. In the formulation developed in the research work presented in this Dissertation, the line models have exactly the same form as the circuits in f i g . 1.3. However, the resistance and equivalent current source are derived in such a way that the frequency dependence of the parameters and the distributed nature of the losses are taken into account. The development of this model is introduced in Chapter 3. In Chapter 2, presented next, some characteristics of time-domain transient simulations are considered. CHAPTER 2 SIMULATION OF ELECTROMAGNETIC TRANSIENTS 2.1 Complete Network Solution: Time Domain vs. Frequency Domain As indicated in Chapter 1, the solution to the line equations can easily be found in the frequency domain, whereas i t s formulation in the time domain requires many simplifications in order to obtain practical results. However, despite the advantages of modelling transmission lines in the fre-quency domain, the solution of a complete system, in which a large variety of conditions and operations are simulated, is more easily formulated in the time domain. As a principle, time-domain formulations have the advantage of deal-ing directly with real physical phenomena instead of with mathematical (some-times "obscure") equations. This philosophical advantage manifests i t s e l f in practice in many forms. For example, one of the major inconveniences in frequency-domain methods is the simulation of sequential operations (time as a variable does not appear explicitly in the frequency domain). An im-portant example of sequential operations is the opening or closing of circuit breakers at specific times, or as a function of the value of certain variables (e.g. when the current f i r s t passes through zero). Another example is the closing or opening of gaps as a function of the magnitude of the voltage. Simulation of switching operations in frequency domain solutions requires going back and forth between frequency and time domains for each change of switch position, thus making the process quite lengthy and cumbersome. On the other hand, simulation of switching operations in the time domain is straightorward. Time domain formulations also greatly f a c i l i t a t e the modelling 19 20 of other non-linear elements, such as surge arresters and saturable reactors and transformers. In these cases, piece-wise linear or non-linear models can easily be implemented. Because of the mentioned factors, most general-purpose programs solve the transients problem directly in the time domain. It is in connection with these programs that the present work has been developed. Dommel's general approach [10] has been used as the basic reference work for time-domain formulations. 2.2 Network Solution in the Time Domain In Dommel's formulation [10], each element in the system is re-presented at each time step t by a resistance in parallel with a current source. The resistances are constant and the values of the current sources are updated according to the conditions at previous time steps. The general form of the network equations at time t is [G][v(t)] = [i(t)] + [ i h ( t ) ] , (2.1) where [G] = constant nodal conductance matrix [v(t)] = vector with (unknown) node voltages at time t [i(t)] = vector with (known) values of external current sources at time t [d.^(t)] = vector with (known from history) values of internal equivalent current sources at time t. The quantities in equation 2.1 are phase quantities. However, as in-dicated in the preceding chapter, the line models are much simpler i f f i r s t derived in terms of decoupled modes. To incorporate the mode models into equation 2.1, the time-mode quantities must be transformed into time-phase quantities. These 21 transformations are considered next. The relationship between phase and mode quantities i s defined in the frequency domain by transformations 1.9 and 1.10: V , = PV (2.2) ph m v ' and V = <>V <2-3> The time-domain form of these relations is obtained by applying the inverse Fourier transform : v , =? - 1{P.V } (2.4) ph m and i , = ? _ 1{Q.I }. (2.5) ph m In the case of completely balanced lines, the line equations can be exactly decoupled by means of real, constant matrices that are independent of the particular line. Examples of these matrices are the Clark Transformation (a, g, 0 components) [3] and the Karrenbauer Transformation [21]. The Karrenbauer Transformation has the advantage that i t can be applied to any number of phases. (This transformation is the one used by Dommel in the "Electromagnetic Transients Program" (EMTP)([24], [38])). The traditional symmetrical components transformation has the computational disadvantage of having complex elements. For transposed lines the above mentioned matrices can decouple the line to a very good degree of accuracy which depends on the type of transposition scheme. In the general case of unbalanced, untransposed lines the modal transformation matrices are a function of the frequency and depend on the 22 particular line configuration. However, i t has been found in studies of this problem (e.g. [22] and [23]) that, especially in the case of single-circuit lines, reasonably accurate results can s t i l l be obtained assuming constant transformation matrices evaluated at some "mid-range" frequency. Assuming that the transformation matrices P and Q are constant, equations 2.4 and 2.5 give V p h ( t > = P Vxn ( t ) ( 2 - 6 ) and ± p h ( t ) " Q ±m ( t ) • <2'7> These equations establish a direct relationship between the time domain quantities of the different line modes and the corresponding phase quan-t i t i e s in the "physical system". A constant diagonal matrix in the time-mode domain (a resistance in the individual mode equivalent circuits) can be transformed into the time-phase domain as follows. With V = R I , m mm from equations 2.2 and 2.3, V = (PRO" 1) I ph m ph and assuming P and Q constant, v , (t) = (PRO" 1) i (t ) , (2.8) ph m ph thus giving, Rph = P RmQ _ 1 ' <2'9> From this relation, in equation 2.1 23 Gph - V" 1 " QRm_1 P _ 1 • (2-10) The vector with the history current sources from the modal equivalent circuits is transformed directly into phase quantities by relation 2.7: i h p h = Q ± h m - t 2 - 1 1 ) CHAPTER 3 THE PROBLEM OF FREQUENCY DEPENDENCE 3.1 General Considerations Two major conclusions can be derived from the general observations made in the preceding chapters: a) The complete line equations, including the frequency variation of the parameters and the distributed nature of the losses, are much more easily formulated in the frequency domain. b) A general transient solution for a complete system is much more conveniently formulated in the time domain. To reconcile both premises, the analytical problem can be stated in terms of developing a formulation, such that, starting from the complete line equations in the frequency domain leads to a time-domain form that can be incorporated into a general time-domain solution. Due to the complexity of the frequency functions involved, the conversion of the line equations from frequency to time domain has to be done numerically. As a result, a series of factors, such as computer time, numerical s t a b i l i t y , and accuracy, become major considerations. As mentioned in the Introduction, the formulations that have been suggested in the past to solve this problem have not yet provided a satisfactory general solution. One of the d i f f i c u l t i e s with some of these formulations is that the elaborate mathematical manipulations required in the process of solution tend to obscure the physical concepts involved in the problem, thus making i t d i f f i c u l t to localize the sources of "unexpected" problems. 24 25 The present work is based on the strong conviction that a formula-tion that facilitates the physical visualization of the concepts involved, such as one based on the concept of an equivalent c i r c u i t , is easier to understand, and thus the possible sources of error easier to assess, than in th case of other "more mathematical" formulations. In this chapter, the evolution ideas associated with equivalent circuit formulations, from the works of Budner (1970, [11]), Snelson (1972, [12]), Meyer and Dommel (1974, [13]), to the new models developed in this work, are presented. In addition to the mentioned works, very useful concepts have been derived from other formula-tions. These formulations w i l l be referred to when these concepts are con-sidered . 3.2 Admittances Formulation: Budner Budner (1974, [11]) presents a very direct circuit-approach to ob-tain a frequency-dependence line model. From the line solution in the fre-quency domain (equations 1.34 and 1.35), V. = cosh(yji) V - Z sinh(y2-)I (3.1) k m e m and I, = i7 sinh(yJi) V - cosh(yii)I , (3.2) K. Z m m c Budner looks at the line as a two-port network ( f i g . 3.1). By nodal analysis the equations of this network are I. = Y. V. + Y. V (3.3) k kk k km m and I = Y , V. + Y V . (3.4) m mk k mm m 26 Fig. 3.1: Single-phase transmission line as a two-port network. From equations 3.1 and 3.2 the admittances are given by Y, . = Y = •=- coth(v£) kk mm Z ' J c and Y = Y km mk csch(Y^) (3.5) (3.6) Taking the inverse Fourier transform of equations 3.3 and 3.4, the following time domain relations are obtained: H(t) = y k k ( t ) * v v(t) + y ^ ( t ) * v„(t) kmv m (3.7) and i m ( t ) = y t n i(t) * v j t ) + y ^ ( t ) * v j t ) 'km kkv m (3.8) In these equations the symbol "*" indicates convolution of the functions involved, (in general, for any two functions f and g, f(t) * g(t) = / f(u).g(t - u)du = / f ( t - u).g(u)du .) (3.9) —OO —00 Due to the form of the convolution integrals, the y-functions in equations 3.7 and 3.8 can be referred to as "weighting functions". After evaluating the convolution integrals in equations 3.7 and 3.8, i (t) and i (t) can be expressed as a function of their corresponding voltages 27 and of values known from previous time steps: V t } • ys v k ( t ) + i k h ( t ) ( 3- 1 0> and i CO = y c v (t) + i (t) , (3.11) m s m —mh ' where y is a constant, and i . , and i . are functions of past values. From s —kh —mh c these relations, an equivalent circuit having the form of the circuit in f i g . 1.3(b) can be directly obtained. One major d i f f i c u l t y in arriving at these equations, however, is the evaluation of the convolution integrals. To calculate the time-domain form of the admittances (weighting functions in the convolution integrals) Budner uses a fast Fourier transform program. The frequency and time domain forms of these admittance functions (for the zero sequence mode of the line studied by Budner) are shown in f i g . 3.2 (reproduced from reference f11]). A number of numerical complications are encountered when inverse-transforming the Y(w) functions. These problems are due in part to the nature of the discrete Fourier transform and, in part, to the form of the Y(OJ) functions. One of the restrictions of inverse Fourier transform methods is that the frequency functions must be band-limited. In order to cope with this problem, Budner limits the upper frequency to a certain value, at the expense of possible loss of detail in the corresponding time domain forms. Apart from this problem, the cyclic nature and the sharp peaks of the Y(w) functions make accurate numerical evaluation of the corresponding y(t) functions f a i r l y d i f f i c u l t . Improved accuracy requires very small AOJ'S and, hence, long computer times. Even assuming that the y(t) functions could be accurately evaluated, the subsequent steps of this formulation s t i l l present other numerical problems. 28 0 80 160 240 320 400 4S0 S60 640 FREQUENCY (HZ) (XIO'l Ca) 5 10 IS 20 25 30 35 40 TIME (MILLISECONDS) Ch) i 10 "T 9 ° 8 x ~ 7 o 6 i 3. 5 w 4 o Z> i h- J o Z » i o 2 I 200 IbO !20 80 _ 40 2 UJ o g -40 S -80 -I 80S! Hsof •200 160 240 320 400 480 56C 640 720 FREQUENCY (HZ) (XI01) (c) (d) Fig. 3.2: Budner's admittance functions. Zero sequence mode (b) y k k ( t ) . (c) Y. fort. (A) v ( d ) (a) Y k k(.) As can be seen from figs. 3.2(b) and (d), the y(t) functions are made up of series of spikes, at f i r s t very sharp and then becoming very f l a t and extending towards in f i n i t e time (these characteristics are even more pronounced for the positive sequence mode, see reference [11]). Because of this form, evalua-tion of the convolution integrals in equations 3.7 and 3.8 presents two important numerical problems: f i r s t , the signals have to be cut-off after some f i n i t e time, and second, a large amount of information can be lost in the convolutions with the very sharp spikes and the very f l a t ones. In general, the convolution process, which has to be performed at each step of the complete network solution, tends to become very time consuming and in-accurate. It is interesting, for the purpose of understanding the physical significance of this formulation, to consider the following particular c o n d i t i o n s . Let i n equations 3.3 and 3.4 v (w) = 0 -* v ( t ) = 0 (end m s h o r t - c i r c u i t e d ) in m and V k ( u j ) = 1 -r v (t) = 6 ( t ) (impulse a p p l i e d at end k) This r e s u l t s i n : Ik<-> = Y k k ( . ) - i k ( t ) = y k k ( t ) and I (co) = Y, (w) ^ km m i m ( t ) - y , m ( t ) ; k n f y k k a n d ykm t h u s r e P r e s e n t currents i k and i i n the system shown i n f i g . 3.3. ik( t)=y kk(t) Ykm(tHm(t) F i g . 3.3: P h y s i c a l i n t e r p r e t a t i o n of Budner's weighting f u n c t i o n s . The form of the y-weighting-functions ( f i g s . 3.2(b) and (d)) can be v i s u a l i z e d from the c o n d i t i o n s i n the system i n f i g . 3.3. The current impulse o r i g i n a t e d at node k needs a c e r t a i n t r a v e l l i n g time to reach node m. Since a time impulse has a uniform frequency d i s t r i b u t i o n , and since d i f f e r e n t frequencies have d i f f e r e n t t r a v e l l i n g times and a t t e n u a t i o n s , when the perturbc t i o n reaches node m the impulse has become a sp i k e . This i s the f i r s t peak of y k m ( f i g . 3.2(d)) (the negative s i g n i s due to the assumed d i r e c t i o n f o r 30 i ) . The wave r e f l e c t e d a t m t r a v e l s back t o k and, a f t e r the c o r r e s p o n d i n g t r a v e l l i n g times and a t t e n u a t i o n s , i t reaches t h i s node, making up the f i r s t s p i k e i n t h e f u n c t i o n y ^ ( f i g . 3.2 ( b ) ) . A f t e r t r a v e l l i n g t w i c e o v e r t h e l i n e , the s p i k e l o o k s now much f l a t t e r : s m a l l e r magnitudes and more sp r e a d o u t . The wave i s r e f l e c t e d a g a i n and t h e same p r o c e s s c o n t i n u e s i n -d e f i n i t e l y . A nother important o b s e r v a t i o n t h a t can be d e r i v e d from the system of f i g . 3.3 i s t h a t the y ( t ) f u n c t i o n s have to be zero f o r t < 0, t h a t i s , they a r e " c a u s a l " time f u n c t i o n s . T h i s f o l l o w s from the f a c t t h a t t h e s e f u n c t i o n s a r e a s s o c i a t e d w i t h r e a l p h y s i c a l v a r i a b l e s , 1 ^ ( 0 and i ^ X m ) , and as such cannot e x i s t b e f o r e the impulse i s a p p l i e d . T h i s c o n s i d e r a t i o n i s v e r y important because i t s e t s a boundary i n the e v a l u a t i o n o f the c o n v o l u t i o n i n t e g r a l s . T h i s i s i l l u s t r a t e d next f o r one of these i n t e g r a l s : CO OO y k k ( t ) * v k ( t ) = v k ( t ~ u ) y k k ( u ) d u = { v k ( t - u ) y k k ( u ) d u • ( 3 ' 1 2 ) I f the lower l i m i t of t h i s i n t e g r a l was l e s s than z e r o , t h i s would imply t h a t v a l u e s o f v ^ f o r times l a r g e r than t would have t o be known,and i t would not be p o s s i b l e t o f i n d t h e r e l a t i o n 3.10 between v ^ ( t ) and i ^ ( t ) . 3.3 T r a v e l l i n g F u n c t i o n s Approach: S n e l s o n ; Meyer and Pommel By a n a l o g y w i t h Bergeron's i n t e r p r e t a t i o n o f the l o s s l e s s , c o n s t a n t -parameter t r a n s m i s s i o n l i n e , S n e l s o n (1972, [12]) d e f i n e s a new s e t o f v a r i a b l e s t o r e l a t e c u r r e n t s and v o l t a g e s a t t h e ends o f the l i n e . T h i s new s e t of v a r i a b l e s l e a d s t o s i m p l e r w e i g h t i n g f u n c t i o n s than those e n c ountered i n Budner's f o r m u l a t i o n . S n e l s o n ' s i d e a i s f u r t h e r d e v e l o p e d by Meyer and Dommel (1974, [13]) and implemented i n t o 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" (EMTP) o f t h e B o n n e v i l l e Power A d m i n i s t r a t i o n (BPA). 31 The new variables are defined as follows: forward travelling functions: f k ( t ) = v k ( t ) + R l i k ( t ) ' (3.13) fJ^ = v J t ) + R i i (t) m m 1 m (3.14) and backward travelling functions: b k(t) = v k(t) - R ; L i k ( t ) , (3.15) b (t) = v (t) - R i (t) n i m 1 m (3.16) For the subsequent development of these equations, R^ could be chosen as any real constant. However, Snelson makes a particular choice and defines R^ as R = lim Z (w) , where Z (to) is the line characteristic impedance (equation 1.33). or*00 c It i s interesting to observe that for a line with constant parameters and no losses ,the forward and backward travelling functions in equations 3.13 to 3.14 are the same as the functions f^ and f^ in d'Alembert's solution (equations 1.42 and 1.43) . Relations 3.13 and 3.14 defining the functions f, and f can be m visualized physically from the system in f i g . 3 .4, where v k(t) = f k ( t ) - R l i k ( t ) and v (t) = f (t) - I L i (t) m m 1 m 32 Fig. 3.4: Equivalent circuit interpretation of the functions f^ and f . The functions and b^ can be interpreted as containing the rest of the information defining "the internal behaviour of the line". To relate the time-domain functions f and b to the exact line solution in the frequency domain, equations 3.13 to 3.16 are transformed into the frequency domain: forward travelling functions: Fk(w) = Vk(co) + R^Coi) , (3.17) F (w) = V (w) + R, I (a>) . (3.18) m m 1 m and backward travelling functions: B k(co) = V k(oj) - R-^Cw) , (3.19) B m ( u ) = V (u) - R I (u) . (3.20) m m 1 m Comparing equations 3.17 to 3.20 with the exact line solution (equations 1.34 and 1.35), i t follows that B. (u) = A.(u) F (u) + A9(a>) F, (co) (3.21) K 1 m 2 k 33 and where B m(u)) = A^oi) F k ( u ) + A2(w) F m ( u ) , (3.22) and A l (.) = ^ (3.23) cosh(y£) + jl— + — sinh(Y£) \ 1 c / A 2(u) = - ^ jsinh( Y£) . A X( W) . (3.24) In equations 3.21 and 3.22 the information regarding the frequency dependence of Z and of Y is carried by the functions A, and A„. This information is c 1 2 transferred into the time domain by transforming equations 3.21 and 3.22 into this domain, thus giving b k(t) = a ^ t ) * f m ( t ) + a 2(t) * f k ( t ) (3.25) and b (t) = a.(t) * f, (t) + a„(t) * f (t) , (3.26) m 1 k z m where a^(t) and a 2(t) are the weighting functions for the time domain convolu-tions . After evaluating the convolutions in equations 3.25 and 3.26, equations 3.15 and 3.16 (representing the "internal behaviour of the line") give at each time step of the solution the following equivalent representa-tions for the line as seen from nodes k and m: V c > = l x \ ( t ) + 4 h ( t > ( 3 - 2 7 ) and = k v J O + (3.28) m m —mh i . , and i . are defined from the past values of the variables. This equivalent —kh —mh r ^ 34 c i r c u i t has the same form as the one i n f i g . 1.3(b), and t h e r e f o r e i t can be d i r e c t l y i n c o r p o r a t e d i n t o Dommel's g e n e r a l t r a n s i e n t s program. E q u a t i o n s 3.27 and 3.28 a r e s i m i l a r t o the ones o b t a i n e d by Budner ( e q u a t i o n s 3.10 and 3.11). They d i f f e r , however, i n the p r o c e d u r e f o r o b t a i n i n g the r e p r e s e n t a t i o n : i n the p r e s e n t case, through the w e i g h t i n g f u n c t i o n s a ^ ( t ) and a 2 ( t ) . As i n Budner's c a s e , the form o f the w e i g h t i n g f u n c t i o n s can be v i s u a l i z e d p h y s i c a l l y by c o n s i d e r i n g some p a r t i c u l a r c o n d i t i o n s . F o r t h i s p urpose, l e t F. (LO) = 1 and F (w) = 0 . (3.29) k m S u b s t i t u t i o n o f these c o n d i t i o n s i n t o e q u a t i o n s 3.21 and 3.22 g i v e s B. ( t o ) = A „ ( t o ) and B (LO) = A. ( t o ) . (3.30) k I m l T r a n s f o r m i n g r e l a t i o n s 3.29 and 3.30 i n t o the time domain, f k ( t ) = 6 ( t ) , f m ( t ) = 0 , (3.31) b, ( t ) = a„(t) and b ( t ) = aAt) . (3.32) k J. m l With c o n d i t i o n s 3.31, the system i n f i g . 3.4 r e s u l t s i n the system shown i n f F i g . 3.5: P h y s i c a l meaning o f the w e i g h t i n g f u n c t i o n s a ^ ( t ) , a 2 ( t ) . f i g . 3.5. S u b s t i t u t i o n o f r e l a t i o n s 3.32 i n t o r e l a t i o n s 3.15 and 3.16 g i v e s 35 a 2(t) = v k ( t ) - R j i j ^ t ) (3.33) and a (t) = v (t) - R i (t) , l m 1 m (3.34) which express a^(t) and a^(t) in terms of the variables in the circuit of f i g . 3.5. It then follows from this circuit that and a 2(t) = 2v k(t) - 6(t) a. (t) = 2v (t) 1 m (3.35) (3.36) In the system of f i g . 3.5, the weighting functions a^(t) and a^(t) are genera-ted by the voltage impulse 6(t) travelling along the line. The resulting form of these functions i s illustrated in f i g . 3.6. (In this figure, the magnitude of the second and subsequent reflexions in a^(t) and a^(t) has been exaggerated for the purpose of ill u s t r a t i o n . In reality, these secondary spikes are very small.) Fig. 3.6: Weighting functions in Snelson, Meyer, and Dommel's formulation. 36 The numerical differences between Budner's and Snelson's weighting functions can be understood from the physical interpretations i n f i g s . 3.3 and 3.5. In Budner's case, the voltage impulse applied at the beginning of the l i n e undergoes f u l l r e f l e c t i o n when i t reaches the s h o r t - c i r c u i t e d r e c e i v i n g end, and s i m i l a r l y when i t comes back to the sending end. In Snelson's case, the l i n e i s terminated at both ends by the impedance R defined as R1 = lim Z (to) . 1 -L to-xo c Therefore, for each frequency component of the applied impulse, the l i n e termination i s much closer to the l i n e c h a r a c t e r i s t i c impedance Z (to) than c a short c i r c u i t i s (the higher the frequency the more t h i s i s t r u e ) . As a r e s u l t , the r e f l e c t i o n s are much smaller and the spikes in the weighting func-tions attenuate much faster i n Snelson's than i n Budner's case. Since the a-functions decay much faster than the y-functions, the corresponding time-domain convolutions at each step of the network solution require much less computer time. Also, the frequency domain form of these functions i s l e s s o s c i l l a t o r y i n Snelson's than i n Budner's case, thus f a c i l i t a t i n g the numer-i c a l evaluation of the corresponding time domain forms. Snelson's formulation of the l i n e equations, as implemented by Meyer and Dommel in the EMTP at BPA, has given s a t i s f a c t o r y r e s u l t s i n many cases of production experience i n the use of t h i s program. Nevertheless, some problems have been encountered. The main source of problems has been the behaviour of the routines at low frequencies (including 60 Hz). In some cases, a f t e r a correct simulation of the i n i t i a l transient period, the routines have had troubles in simulating the subsequent stationary conditions. Experience in t r y i n g to correct these problems seems to indicate that they are associated with the d i f f i c u l t y i n evaluating with accuracy the " t a i l " portion of the a^(t).and a^(t) functions (the portion of these functions a f t e r the main spikes). Accurate evaluation of these t a i l portions i s important because despite th e i r 37 small magnitudes the spikes extend over a long period of time,and therefore their total contribution to the convolution integrals can be significant. In Meyer and Dommel's implementation, the time domain form of the weighting functions is calculated from their corresponding frequency domain definitions (equations 3.23 and 3.24) through inverse Fourier transforms. Comparisons in which the time domain functions are brought back into the frequency domain by the direct Fourier transform have given smaller magni-tudes than the original values in the low frequency region. These results seem to corroborate the erratic behaviour of the routines at low frequencies. Also, in these studies, the function a 2(t) has presented more error than the function a^(t). 3.4 New Formulation Several considerations lead to the belief that, as mentioned at the end of the preceding section, the d i f f i c u l t i e s in evaluating with accuracy the t a i l portions of the a-weighting-functions in Meyer and Dommel's formula-tion are responsible for the problems encountered with this method at low frequencies. The resistance of the line (especially for the zero sequence mode) increases considerably with frequency. This means that the higher frequency components of the signals are attenuated much faster than the lower frequency components. Therefore, in the t a i l portions of a^(t) and a 2(t) there is much less contribution of high frequencies than of low frequencies. In addition to this factor, the resistance in f i g . 3.5 is much closer to the line characteristic impedance at high frequencies than i t is at low frequencies. This means that the more reflections the applied voltage impulse undergoes, the higher the ratio of low frequency components with respect to high frequency components. 38 A c c o r d i n g to the p r e c e d i n g c o n s i d e r a t i o n s , one way of i m p r o v i n g the low f r e q u e n c y r e s p o n s e would be t o choose f o r a v a l u e t h a t i s c l o s e r t o the c h a r a c t e r i s t i c impedance f o r low f r e q u e n c i e s . But, t h i s would p r o b a b l y worsen the h i g h f r e q u e n c y r e s p o n s e and, c o n s e q u e n t l y , the s i m u l a t i o n o f t h e v e r y f a s t i n i t i a l p e r i o d o f the t r a n s i e n t phenomenon. A compromise would be to make R^ c l o s e to Z^Cco) f o r some mid-range f r e q u e n c y . F o l l o w i n g the above l i n e o f r e a s o n i n g , a much more d e f i n i t i v e s o l u t i o n t o t h e problem would be t o have a t the ends o f the l i n e i n f i g . 3.5, i n s t e a d o f a s i m p l e r e s i s t a n c e R^, a c e r t a i n network whose f r e q u e n c y response matches the l i n e c h a r a c t e r i s t i c impedance, Z^(us) , over the e n t i r e f r e q u e n c y range. T h i s i s shown i n f i g . 3.7. In t h i s case the v o l t a g e impulse g e n e r a t e d 5(t) Network Z eq ik(t) line vk(t) a2(t) = o im«) a,(t) ——rrt •vm(t)/ 1 Network Zeq F i g . 3.7: P h y s i c a l i n t e r p r e t a t i o n o f the new w e i g h t i n g f u n c t i o n s . a t node k would not be r e f l e c t e d back from node m, and a ^ ( t ) would c o n s i s t of o n l y one s p i k e ; a 2 ( t ) would be c o m p l e t e l y e l i m i n a t e d , t h a t i s , a 2 ( t ) = 0. The r e s u l t i n g new w e i g h t i n g f u n c t i o n s would thus have the form shown i n f i g . 3.8. How to f o r m u l a t e the l i n e e q u a t i o n s i n o r d e r to a r r i v e a t t h e w e i g h t i n g f u n c t i o n s g e n e r a t e d by the system i n f i g . 3.7 i s c o n s i d e r e d n e x t . 39 (a) eg o zero (b) Fig. 3.8: Weighting functions in the new formulation. (a) a 1 ( t ) . (b) a 2(t) (equal to zero). General equations: For the new formulation, Snelson's travelling functions are redefined as follows: forward travelling functions: f k ( t ) = v k(t) + e k(t) , f (t) = v (t) + e (t) , m m m (3.37) (3.38) and backward travelling functions: b k(t) = v k(t) - ek(.t) , b (t) = v ( t ) - e (t) m m m (3.39) (3.40) 40 where e, (t) = voltage across equivalent network Z due to i,(t) k eq K and e (t) = voltage across equivalent network Z due to i (t) . m b ^ eq m Equations 3.37 and 3.38 are represented in the system shown in f i g . 3.9. Equations 3.39 and 3.40 contain the information leading to the new weighting functions. ik(t) fk(t)@ I l i n e VJt) N etwork Z e q ek(t) ± ) W « Network Z e q Fig. 3.9: Equivalent circuit interpretation of the new f, and f functions, k m Assuming that the network Z is linear, equations 3.37 to 3.40 can be transformed into the frequency domain as follows: forward travelling functions: Fk(w) = V k(uj) + Z (o>) I k ( o i ) , (3.41) F (oo) = V (to) + Z (to) I (to) , (3.42) m m eq m v ' and backward travelling functions: s B k ( u ) = V k(to) - Z (u) I k ( u ) , (3.43) B m ( « ) = V (cu) - Z o (u) I (to) . (3.44) m m eq m These equations correspond to Snelson's 3.17 to 3.20 substituting by Z (to) The voltages across Z in the time domain are related to the corresponding voltages in the frequency domain by e f c ( t ) « • E k(to) = Z (u>) I k ( u > ) (3.45) and if e m ( t ) < • E m(o J) = Z e q ( . ) I m ( . ) , (3.46) where Z (to) = Z (to) (3.47) eq c (within the limits of the approximation). Comparing equations 3.41 to 3.44 with the general solution for the line equations in the frequency domain, equations 1.34 and 1.35, and with 3.47, i t follows that B. (to) = A, (w) F (to) (3.48) k i m and B m(to) = A 1(to) F k(co) , (3.49) where - Y ( C O ) I 1 A 1 ( W ) = 6 = cosh( Y*) + B l n h ( Y * ) ' ( 3 - 5 0 ) Referring back to Meyer and Dommel's A^(to) and A^iu) functions (equations 3.23 and 3.24), i t can be seen that by replacing in these equations R^ by Z c(co) , equation 3.50 is obtained for A^(co) and A^iw) = 0, thus corresponding to the results of the present derivation. The time domain form of equations 3.48 and 3.49 is given by 42 b k ( t ) = a x ( t ) * f m ( t ) (3.51) and b m ( t ) = a x ( t ) * f k ( t ) , (3.52) which correspond to equations 3.25 and 3.26 i n Meyer and Dommel's formulation. The convolution i n t e g r a l s i n equations 3.51 and 3.52 are given by b, (t) = / f (t - u) a n(u)du = r f (t - u) 3 l ( u ) d u (3.53) K. m l m l and b m ( t ) = f" f f c ( t - u) a ] L(u)du = / f (t - u) a (u)du . (3.54) —OO The lower l i m i t i n the above i n t e g r a l s i s x, since a^(t) = 0 for t < x ( f i g . 3.8). With reference to f i g . 3.7, x represents the time needed by the fas t e s t frequency component of the applied voltage impulse to reach the receiving end of the l i n e . (This time i s given by x = l i n e length . /LC for to •+ <*>.) In a d i s c r e t e numerical s o l u t i o n of the system network, i f the time step At i s chosen to be smaller than x, then the values for f, and f i n the k m convolution i n t e g r a l s 3.53 and 3.54 are known (at d i s c r e t e points) from the previous time steps and these i n t e g r a l s can be evaluated. Since the new weighting function a^(t) consists of only one spike, the upper l i m i t i n the "history i n t e g r a l s " , equations 3.53 and 3.54, does not extend beyond x for very long, and the time i n t e r v a l during which these i n t e -grals have to be evaluated i s much shorter than i n the case of Meyer and Dommel's formulation. From the preceding considerations, the following advantages of the new formulation compared to Meyer and Dommel's can be noted: a) The problem of accurate evaluation of the t a i l portions of the weighting functions has been eliminated. 43 b) The i n t e r v a l of i n t e g r a t i o n to e v a l u a t e the h i s t o r y f u n c t i o n s i s much s h o r t e r , thus c o n s i d e r a b l y r e d u c i n g computer time i n n u m e r i c a l e v a l u -a t i o n s of these i n t e g r a l s . c) The problem of a c c u r a t e e v a l u a t i o n of the w e i g h t i n g f u n c t i o n a^(t) has been e l i m i n a t e d . P h y s i c a l i n t e r p r e t a t i o n o f a]_(t) i n the new f o r m u l a t i o n : The o r i g i n a l o b j e c t i v e i n e s t a b l i s h i n g the new f o r m u l a t i o n was to a r r i v e at the p h y s i c a l i n t e r p r e t a t i o n i n f i g . 3.7. That t h i s o b j e c t i v e has been accomplished i s v e r i f i e d n e x t . L e t i n e q u a t i o n s 3.48 and 3.49 F ( t o ) = 0 and F. (OJ) = 1 . (3.55) m K. Then B ( t o ) = 0 and B (OJ) = A, ( t o ) , (3.56) k m l which i n the time domain g i v e s f m ( t ) = 0 , f k ( t ) = 6 ( t ) , (3.57) b. ( t ) = 0 and b ( t ) = a , ( t ) . (3.58) k m i With c o n d i t i o n s 3.57 i n the o r i g i n a l d e f i n i t i o n of the new forward t r a v e l l i n g f u n c t i o n s ( e q u a t i o n s 3.37 and 3.38), i t f o l l o w s t h a t v k ( t ) = 6 ( t ) - e k ( t ) (3.59) and v ( t ) = - e ( t ) . (3.60) m m With r e l a t i o n s 3.57, 3.59, and 3.60, the system of f i g . 3.9 becomes the o r i g i n a l l y d e s i r e d system o f f i g . 3.7. S u b s t i t u t i o n of c o n d i t i o n s 3.58 i n 44 the new backward t r a v e l l i n g f u n c t i o n s ( e q u a t i o n s 3.39 and 3.40) g i v e s V k ( t ) = e k ( t ) ( 3 - 6 1 ) and v m ( t ) = a i ( t ) + e m ( t ) . (3.62) From e q u a t i o n s 3.60 and 3.62 i t then f o l l o w s t h a t a l ( t ) = 2 v m ( t ) ' (3.63) thus verifying the original assumption regarding the form of the new a ^ t ) weighting function. The function a„(t) does not appear in this new formulation. L i n e e q u i v a l e n t c i r c u i t : At a g i v e n time s t e p t o f the network s o l u t i o n , the v a l u e s o f the f u n c t i o n s b. ( t ) and b ( t ) a r e d e f i n e d from the h i s t o r y o f the l i n e v a r i a b l e s k m J ( e q u a t i o n s 3.53 and 3.54). S u b s t i t u t i o n o f these v a l u e s i n the o r i g i n a l d e f i n i t i o n o f the f u n c t i o n s b ^ and b ^ ( e q u a t i o n s 3.39 and 3.40) g i v e s the l i n e e q u i v a l e n t c i r c u i t s a t nodes k and m. That i s , l e t , from e q u a t i o n s 3.53 and 3.54, b ^ ( t ) = e ^ (from h i s t o r y ) (3.64) and b ( t ) = e , (from h i s t o r y ) . (3.65) m —mh Then, from e q u a t i o n s 3.39 and 3.40, v k ( t ) = e k ( t ) + e^Ct) (3.66) and v ( t ) = e (t) + e , ( t ) , (3.67) m m —mh 45 where e and e represent the voltages across the network Z . Equations 3.66 K m eq and 3.67 can be represented by the equivalent circuits shown in f i g . 3.10. N e t w o r k i_K(t) Z, eq £ Z Z > + ek(t)-v„(t) N e t w o r k Zeq im(t) m © e k h ( t ) fimh(t) em(t) + f vm(0 (a) 00 Fig. 3.10: New line models at time step t. (a) Node k. ( e ^ and known from past values .) (b) Node m. The circuits in f i g . 3.10 have the same form as the cir c u i t in f i g . 1.3(a) of Dommel's representation for the simplified line. The only difference in the circuits is that the simple resistance R , representing the character-i s t i c impedance of the simplified line, has now been replaced by the network Z^, which represents the frequency-dependent characteristic impedance of the real line. As w i l l be discussed in detail later in this work, the network Z eq in the models in f i g . 3.10 can be simulated, to a high degree of accuracy, by a combination of Resistance-Capacitance (R-C) building blocks. As in Dommel's general transient program [10], these R-C blocks can in turn be simulated by pure resistances in combination with corresponding history current sources. The line models in f i g . 3.10 can then be further reduced to an equivalent resistance in parallel with a total equivalent history current source, and thus to the same simple form of the circuit in f i g . 1.3(b). The idea of using R-C combinations to model the line characteristic impedance is considered in the work by Groschupf [25], who uses a simple 46 R-C model to approximate Z^Cco) over the low-and mid-frequency ranges. Yet, Groschupf concludes that the effect of the frequency dependence of the characteristic impedance is not very significant. The results presented later in this Dissertation, however, show that Groschupf's conclusions are correct only for some simple cases of open-circuited lines. For a general class of transient studies, and particularly when short circuits are involved, the correct simulation of Z^Cw) can become very important. Some additional physical considerations: In order to explore further the physical meaning of the weighting function a^(t), i t is interesting to consider the simplified case when the system has no losses and the parameters do not depend on the frequency. Under these conditions there w i l l be no distortion or attenuation, and in the system in f i g . 3.7 the impulse injected at node k w i l l arrive unmodified at node m after a time x. The comparison between the form of a^(t) for this ideal case and the form of a.,(t) for the general case is shown in f i g . 3.11. (a) (b) Fig. 3.11. Weighting function a^(t): (a) general case; (b) ideal case. For the ideal a^(t) in f i g . 3.11(b), the history function D m ( t ) in equation 3.54 becomes 47 CD b m ( t ) = / f k ( t - u) a ] L ( u ) d u = f k ( t - u) 6(u - T ) d u = f k ( t - T ) , (3.68) and w i t h e q u a t i o n s 3.40 and 3.37, v m ( t ) - e m ( t ) = v k ( t - T ) + e k ( t - T ) . (3.69) But, f o r c o n s t a n t parameters and no l o s s e s Z (OJ) = R = c o n s t a n t , c c and then e ( t ) = R i ( t ) (3.70) • m e m and e k ( t - x) = R c i k ( t - T ) , (3.71) f i n a l l y g i v i n g , v ( t ) - R i ( t ) = v, ( t - x) + R i , ( t - x) . (3.72) m c m k c K. T h i s e q u a t i o n i s e x a c t l y the same as e q u a t i o n 1.46 o f Dommel's f o r m u l a t i o n f o r the s i m p l i f i e d l i n e r e p r e s e n t a t i o n . Comparing the forms o f a ^ ( t ) f o r the g e n e r a l and s i m p l i f i e d cases ( f i g . 3.11), i t i s i n t e r e s t i n g to no t e t h a t i n t h e g e n e r a l case t h e r e l a t i o n -s h i p between c u r r e n t and v o l t a g e at one o f the nodes a t time t depends not o n l y on the p a s t c o n d i t i o n s at the o t h e r node x u n i t s o f time e a r l i e r (as i n the i d e a l c a s e ) , b u t on the "weighted" e f f e c t o f p a s t c o n d i t i o n s from x to the time at which a ^ ( t ) becomes n e g l i g i b l e . T h i s i s a consequence o f d i f f e r e n t f r e q u e n c y components h a v i n g d i f f e r e n t t r a v e l l i n g times and a t t e n u a t i o n s . The p r e c e d i n g o b s e r v a t i o n s i l l u s t r a t e the e f f e c t o f a " n o n - i d e a l " w e i g h t i n g f u n c t i o n . The e f f e c t o f a " n o n - i d e a l " c h a r a c t e r i s t i c impedance ( t h a t i s , o f a frequency-dependent Z^) can be seen from e q u a t i o n 3.69, where a ^ ( t ) has been assumed to be i d e a l , b u t t h e v o l t a g e s a c r o s s the e q u i v a l e n t 48 impedance, e^ and e^, are s t i l l g e n e r a l . F o r the e q u i v a l e n t h i s t o r y v o l t a g e s o u r c e , e q u a t i o n 3.69 g i v e s the same model as Dommel's s i m p l i f i e d r e p r e s e n t a t i o n i n f i g . 1.3(a). The o n l y d i f f e r e n c e i s t h a t now the s i m p l e r e s i s t a n c e R i s c r e p l a c e d by the e q u i v a l e n t c h a r a c t e r i s t i c impedance network S i z e o f the i n t e g r a t i o n s t e p : F o r the e v a l u a t i o n of the h i s t o r y c o n v o l u t i o n i n t e g r a l s i n e q u a t i o n s 3.53 and 3.54, i t was assumed t h a t the i n t e g r a t i o n s t e p At i n the complete network s o l u t i o n was s m a l l e r than the l i n e t r a v e l l i n g time x. However, i n some s t u d i e s combining r e l a t i v e l y l o n g and r e l a t i v e l y s h o r t l i n e s , i t might seem c o n v e n i e n t , i n o r d e r to save computer t i m e , t o choose a At t h a t i s l a r g e r than the t r a v e l l i n g time of some of the s h o r t e r l i n e s . A p o s s i b l e s o l u t i o n to t h i s problem would be to r e p r e s e n t the s h o r t e r l i n e s by models t h a t i n c l u d e the e f f e c t of x < At. T h i s case i s mentioned i n Meyer and Dommel's work of r e f e r e n c e [13] and l e a d s to mutual c o u p l i n g between the e q u i v a l e n t c i r c u i t s a t nodes k and m. T h i s c o u p l i n g c o m p l i c a t e s the s o l u t i o n of the complete system network. Another a l t e r n a t i v e to the m o d e l l i n g of " r e l a t i v e l y s h o r t " l i n e s would be to r e p r e s e n t these l i n e s by frequency-dependent lumped-parameter s e c t i o n s . The m o d e l l i n g of f r e q u e n c y dependence i n lumped-parameter elements i s b r i e f l y d i s c u s s e d n e x t . Frequency dependence of lumped-parameter elements: In t r a n s i e n t s t u d i e s many of the system components a r e r e p r e s e n t e d by lumped-parameter models. F o r many o f these components (e . g . power t r a n s -f o r m e r s , r e a c t o r s , g e n e r a t o r s , and l o a d s ) a more a c c u r a t e r e p r e s e n t a t i o n c o u l d be o b t a i n e d by t a k i n g i n t o account the f r e q u e n c y dependence o f t h e i r parameters. With the techniques developed in this work for the modelling of frequency-dependent transmission lines, the modelling of frequency-dependent lumped-parameter elements can be considered as a particular, simplified case. In the equivalent circuits of f i g . 3.10 there are two phenomena &o be modelled: the network Z representing the frequency-dependent character-i s t i c impedance, and the history voltage sources representing the voltages and currents travelling along the line from the opposite end. In the case :of lumped parameters there are no travelling times involved, and the problem i<s reduced to the modelling of a frequency-dependent impedance. The techniques described later in this work for the modelling of the line characteristic impedance can be directly applied, as well, to the modelling of the frequency-dependent impedance of a lumped element. The application of these concepts for specific cases of frequency-dependent lumped-parameter elements is beyond the scope of this Dissertation and is le f t as an area of future research. CHAPTER 4 SYNTHESIS OF THE CHARACTERISTIC IMPEDANCE AND WEIGHTING FUNCTION 4.1 G e n e r a l C o n s i d e r a t i o n s The p r o c e d u r e s f o r o b t a i n i n g the parameters of the new f r e q u e n c y -dependence l i n e model are d i s c u s s e d i n t h i s c h a p t e r . S i n c e the e q u i v a l e n t c i r c u i t s a t b o t h l i n e ends are analogous t o each o t h e r , i n o r d e r to s i m p l i f y the d e s c r i p t i o n , o n l y one of the ends i s c o n s i d e r e d . E q u a t i o n s 3.43 and 3.39 g i v e the b a s i c r e l a t i o n s h i p s f o r the l i n e e q u i v a l e n t c i r c u i t at node k i n the f r e q u e n c y and time domains: B k ( W ) = Vk(o>) - Z e q(co) I k ( u i ) (4.1) and b k ( t ) = v k ( t ) - e k ( t ) . (4.2) The c o r r e s p o n d i n g e q u i v a l e n t c i r c u i t s a r e shown i n f i g . 4.1. I k ^ Z e q M 1 f © Bk(«J Network ' k< 0 Z e q k P H t-Vic(t) ® b k ( t ) (a) (b) F i g . 4.1: L i n e e q u i v a l e n t c i r c u i t a t node k. (a) Frequency domain. (b) Time domain. 50 51 i In f i g . 4.1(b) the network Z is a constant-parameter linear net-eq work whose frequency response Z (OJ) simulates the line characteristic im-eq pedance, that i s , ideally, Z (ID) = Z (w). The history voltage source b, (t) ecj c K, is given by equation 3.53 as CO b, (t) = / f (t - u) a i(u)du. (4.3) K. x m 1 In this integral, a^(u) is the history weighting function and is given by the time domain form of equation 3.50. The form of a^(t) was illustrated in f i g . 3.8(a). The function f (t - u) represents the past values of the voltage m and current at the other end of the line (node m in this case) and is given, from i t s definition in equation 3.38, by f (t - u) = v (t - u) + e <t - u) . (4.4) m m m In integral 4.3 the most recent value of f corresponds to u = T , that i s , f (t - x). As u increases, and for as long as a,(u) ^ 0, the values of f m 1 m go further back in time. For u > t, the values of f correspond to times before m the transient process was started, that i s , to i n i t i a l conditions. In equation 4.4 the values of v^ are the actual voltages at node m, and are thus known from the solution of the network at previous time steps and from i n i t i a l conditions. The values of e can be obtained from the past m values of v and from the past values of relation 3.40: m e (t - u) = v (t - u) - b (t - u) . (4.5) m m m Actually, e (t - u) from this relation can be substituted in equation 4.4, m giving a more direct relation for f (t - u): f (t - u) = 2v (t - u) - b (t - u) . (4.6) m m m 52 The e v a l u a t i o n of the parameters of the c i r c u i t i n f i g . 4.1(b) p r e s e n t s , t h u s , two b a s i c a s p e c t s : a) S y n t h e s i s o f a l i n e a r network to s i m u l a t e the l i n e c h a r a c t e r -i s t i c impedance. b) C a l c u l a t i o n of the w e i g h t i n g f u n c t i o n a ^ ( t ) and e v a l u a t i o n of the h i s t o r y v o l t a g e source b ( t ) . K. These problems are d i s c u s s e d n e x t . 4.2 S y n t h e s i s of E q u i v a l e n t Network to Represent the C h a r a c t e r i s t i c Impedance 4.2.1 G e n e r a l C o n s i d e r a t i o n s From e q u a t i o n 1.33 the c h a r a c t e r i s t i c impedance f o r a g i v e n l i n e mode i s g i v e n by „ , v /Z' (OJ) - , N Vw) = /rfcy> (4-7) where, i n g e n e r a l , Z'(w) = R'(u) + ja)L ' ( u ) ) (4.8) and Y'(w) = G'(o)) + ju)C'(w) . (4.9) The v a r i a t i o n o f the s e r i e s r e s i s t a n c e (R') and i n d u c t a n c e (L') w i t h f r e q u e n c y i s i l l u s t r a t e d i n f i g . 4.2. The c u r v e s i n t h i s f i g u r e c o r r e s p o n d to the z e r o and p o s i t i v e sequence modes of a t y p i c a l 100-mi, 500kV, 3-phase t r a n s m i s s i o n l i n e , and were c a l c u l a t e d u s i n g the UBC L i n e C o n s t a n t s Program [ 2 6 ] . T h i s program uses Carson's formulae [19] to e v a l -u a t e the f r e q u e n c y dependence of the l i n e parameters. 53 io- q IO1 = 10' d io-; 10- i II II 11 11 IM n—11 MI M—i min—i mm—i mm—11 in n—i null—11 in u—11 HI n—i nun 10"' 10 10" 1 10 10' 103 10* 10s 10* JO1 10* FREQUENCY (HZ) 0 -\—11 HI II i mi II—i i ill ii—i IIIIII—i 11 nil—i i in n—r io io io i io io! io io* )if io io itr FREOUENCT (HZ) (a) (b) F i g . 4.2: L i n e parameters R' and L' as a f u n c t i o n o f f r e q u e n c y (a) R' z e r o and p o s i t i v e sequence. (b) L' z e r o and p o s i t i v e sequence. In a c t u a l t r a n s m i s s i o n l i n e s the c a p a c i t a n c e C' i s not v e r y s e n s i t i v e to v a r i a t i o n s i n f r e q u e n c y , and f o r f r e q u e n c i e s up t o about 1 MHz (e.g. Hedman [9]) i t i s u s u a l l y assumed to be c o n s t a n t . The shunt conductance G' i s a s s o c i a t e d m a i n l y w i t h l e a k a g e a t the i n s u l a t o r s , and i t s v a l u e i s a f f e c t e d by a s e r i e s o f e x t e r n a l f a c t o r s , such as h u m i d i t y , p o l l u t i o n , a n d maintenance. The c o r r e c t v a l u e o f G' i s t h e r e f o r e d i f f i c u l t t o a s s e s s w i t h a c c u r a c y . F o r t u n a t e l y , t h i s v a l u e i s n o r m a l l y s m a l l compared t o the o t h e r system parameters ( e . g . i n the o r d e r o f .3 x 10 7 ^/km, Dommel [27] p. 57), and i t s e f f e c t upon t h e system t r a n s i e n t c o n d i t i o n s i s n o r m a l l y not v e r y s i g n i f i c a n t . However, as e x p l a i n e d n e x t , i t i s c o n v e n i e n t f o r numer-i c a l r e a s o n s to c o n s i d e r a f i n i t e v a l u e o f G' i n the development o f the l i n e t r a n s i e n t e q u i v a l e n t c i r c u i t s . From e q u a t i o n 4.7, f„\ /R'(o)) + JML'(tt)) >> - / G« + jcoC (4.10) 54 Fo r 03 = 0 (dc c o n d i t i o n s ) , s i n c e R'(0) and L'(0) a r e f i n i t e , Z (dc) c (dc) (dc) (4.11) I f G'(dc) i s assumed t o be z e r o , then Z c ( d c ) ->- 0 0. S i n c e Z c ( to ) appears as a f a c t o r ( m u l t i p l y i n g o r d i v i d i n g ) i n the f o r m u l a t i o n o f the l i n e e q u a t i o n s , «"> l.| HI" ) io io io 10 10* I01 10' 10" Iff 10' icf FREQUENCY IHZI FRE0UFNCT IHZI (a) (b) .... . nimi i nun i HUM I nmi i mm icr icr io itr icr 10 itr FREDUENCT IHZI -100 1 limn iiimi limn limn nimi nimi nnin Minn iium nimi I nim io io io" i io irf icr io" itf irf 10 itr rREOUENCT IHZI (c) (d) F i g . 4.3: C h a r a c t e r i s t i c impedance Z c ( t o ) . Zero sequence mode. a) Magnitude, b) Phase a n g l e , c) R e a l p a r t , d) Imaginary p a r t , 55 having Zc(io) •*• °° for small or zero frequencies can cause numerical i n s t a b i l i t y problems in transient studies involving d.c. or very low frequency components. The form of the functions Z (u) for the zero and positive modes of the c reference line mentioned before are shown in figs. 4.3 and 4.4. These functions were evaluated from equation 4.10, where R1, L', and C'were obtained from the UBC Line Constants Program and G'was assumed to be equal to .3 x 10 7 /^km. io io io i io io" \€ io iof itr to irr rREOUENCt CHZI 30 25 20 _ 15 8 io ? u ~ -10 -J £ -.5 -20 -25 -30 10 to" io" 10* 101 10* 10' 10* FREQUENCY IMZ1 (a) (b) 200 I inmi iniui I I I H I i ni.u i nun i min 11 iimii mini iinni to io io i io iff i<r io irf iir io i<r FREQUENCY IHZI 200 150 100 50-1 i _ 5 H : .ioo -150 -200 I 11inn 11IBn i nun Ii 10'" 10* 10'' I 10 IO1 10" 10" 10* 10* 10' 10* FREOUENCT IHZI (c) (d) Fig. 4.4: Characteristic impedance Z c (cu) . Positive sequence mode. a) Magnitude, b) Phase angle, c) Real part, d) Imaginary part. 56 4.2.2 Rational Approximation of Z C(CJ) The problem of finding a linear network whose frequency response approximates Z^COJ) is a problem of network synthesis theory. A series of classical methods exist in network synthesis theory (e.g. Guillemin [28]) for the realization, using R, L, and C constant elements, of a given rational impedance function (a function expressed as a fraction of polynomials in the complex variable s, s = a + ju>) . In the present case, however, the function Z c(oo) is not known as a rational function, but only as a point by point (tabular) function of frequency. Therefore, the f i r s t step in order to simu-late Z c(w) by an equivalent linear network i s to find a rational approximation for this function. Because of physical and numerical considerations that w i l l be explained later, the rational function to simulate Z^ito) is chosen to have the following general form: (s + z,)(s + z„) ... (s + z ) Z (s) = H . . \ , ^ Z . , . " , (4.12) eq v (s + p 1)(s + p 2) ... (s + p ) ' with the following characteristics: (a) The constant H is real and positive. (b) The number of poles i s equal to the number of zeros. (c) The poles and zeros are real, negative, and simple. It can directly be seen that with this definition of Z the boundary condi-J eq J tions of Z (to) (figs. 4.3 and 4.4) can be satisfied. With s = jto and i) OJ ->• 0 Z 1 Z 2 * " Zn Z (s) -»- H = real positive constant, eq P j P 2 ••• P n which should correspond to Z (dc) in equation 4.3. c 57 i i ) ze q ( s ) H = r e a l p o s i t i v e c o n s t a n t which s h o u l d c o r r e s p o n d to Z c ( i o ) i n e q u a t i o n 4.10 f o r to 4.2.3 E q u i v a l e n t R-C Network A r a t i o n a l f u n c t i o n w i t h r e a l , n e g a t i v e , and s i m p l e p o l e s and z e r o s and w i t h the number of p o l e s and z e r o s not d i f f e r i n g by more than one can be r e a l i z e d by an R-C ( r e s i s t a n c e - c a p a c i t a n c e ) e q u i v a l e n t network. The f u n c t i o n Z e q ( s ) i n e q u a t i o n 4.12 s a t i s f i e s these c o n d i t i o n s . S i n c e the l i n e e q u i v a l e n t model i n f i g . 4.1(b) i s f i r s t f o r m u l a t e d i n s e r i e s form, i t i s c o n v e n i e n t f o r the f u r t h e r r e d u c t i o n of t h i s c i r c u i t t o s y n t h e s i z e Z by a s e r i e s network. eq F o s t e r I r e a l i z a t i o n : The F o s t e r I , R-C s y n t h e s i s ( e . g . K a m i [ 2 9 ] , p. 154) i s based on a s i m p l e p a r t i a l f r a c t i o n e x p a n s i o n of the network f u n c t i o n . In the case of e q u a t i o n 4.12, Z^(s) can be e x p r e s s e d as k, k 0 k Z (s) = kr, + — + - — + + n e q 0 s + P l T s + p 2 + + s + p » n (4.13) where and k 0 = l i m Z e q ( s ) = H ' (4.14) k = (s + p ) Z (s) i x eq s = - p, (4.15) ( i = 1, 2, n) ( I n t h i s l a s t r e l a t i o n , m u l t i p l i c a t i o n by the f a c t o r (s + p.) i s to be u n d e r s t o o d as a c a n c e l l i n g out o f the c o r r e s p o n d i n g f a c t o r i n the denominator 58 of Z (s) , or as the l i m i t of the product when s -> - p..) eq x In equation 4.13 the f i r s t term of Z e ^ ( s ) corresponds d i r e c t l y to a r e s i s t a n c e R, and the other terms correspond to R-C p a r a l l e l combinations: R // -7; 1/C sC s + 1/RC (4.16) The network corresponding to equation 4.13 i s , t h e r e f o r e , as shown i n f i g . 4.5. The parameters of t h i s network are given by: and R 0 " k o (4.17) R. = k./p. C. = 1/k. (4.18) (4.19) f o r i = 1, 2, , n, R 2 F i g . 4.5: Foster I r e a l i z a t i o n of Z eq I t can r e a d i l y be seen that the equivalent network i n f i g . 4.5 can f o l l o w the behaviour of the l i n e c h a r a c t e r i s t i c impedance ( f i g s . 4.3 and 4.4) at the boundary c o n d i t i o n s . For very low frequencies the capacitances tend to open c i r c u i t s , and at the l i m i t (OJ = 0 (dc c o n d i t i o n ) ) Z (o> = 0) = Z (dc) = R n + R + R„ + ... + R eq c 0 1 2 n (4.20) For high frequencies the capacitances tend to short c i r c u i t s , and f o r OJ •> 0 0 Z (w °°) = Z (w -»• o°) = R„ eq c 0 (4.21) 59 4.3 W e i g h t i n g F u n c t i o n and H i s t o r y C o n v o l u t i o n 4.3.1 G e n e r a l C o n s i d e r a t i o n s A f t e r f i n d i n g an e q u i v a l e n t network to s i m u l a t e the c h a r a c t e r i s t i c impedance, the next s t e p i n c o m p l e t i n g the l i n e e q u i v a l e n t c i r c u i t o f f i g . 4.1(b) i s to f i n d t h e h i s t o r y v o l t a g e s o u r c e , b ^ ( t ) . T h i s f u n c t i o n i s g i v e n by 00 b : ( t ) = / f ( t - u) a 1 ( u ) d u . (4.22) k m 1 To e v a l u a t e t h i s i n t e g r a l , the w e i g h t i n g f u n c t i o n a ^ ( t ) has f i r s t to be determined. The f r e q u e n c y domain form o f t h i s f u n c t i o n , A^(to), i s known from i t s d e f i n i t i o n i n terms of the system parameters. From e q u a t i o n 3.50, A 1(co) = e , (4.23) where £ i s the l i n e l e n g t h , and y ( a ) ) 1 S g i v e n from e q u a t i o n 1.32 by Y(UJ) =/Z'(W) Y'(co) , (4.24) which, w i t h the c o n s i d e r a t i o n s made e a r l i e r i n c o n n e c t i o n w i t h the f r e q u e n c y dependence o f the parameters, can be e x p r e s s e d as Y(UJ) = /[R'(UJ) + j o j L ' M n G * + jtoC'] . (4.25) The form o f A^(co) f o r the ze r o and p o s i t i v e sequence modes of the r e f e r e n c e 100-mi l i n e i s shown i n f i g s . 4.6 and 4.7. With A^(oj) e v a l u a t e d from the system p a r a m e t e r s , the f u n c t i o n a ^ ( t ) can be de t e r m i n e d by means of an i n v e r s e F o u r i e r t r a n s f o r m n u m e r i c a l t e c h n i q u e . Once a ^ ( t ) has been determined f o r a p a r t i c u l a r l i n e mode, the h i s t o r y f u n c t i o n i n e q u a t i o n 4.22 can a l s o . b e e v a l u a t e d a t each time s t e p of the network s o l u t i o n by a n u m e r i c a l i n t e g r a t i o n t e c h n i q u e . There i s an 60 io io1 ity io* \tf io" io' FREQUENCY IHZI 180 • ISO • 120-90 • 8 GO • ai IU «-> 30 • i 0 --30 • a -60 • -90 • -120 • -ISO • -160 • 10'' I 10 10" FREQUENCY IHZI 10* (a) (b) ID*1 2 4 1 2 4 10 2 4 IO* 2 4 10* Z FREQUENCY IHZI 10" 2 4 ID* 2 4 10 2 4 10" 2 4 10* 2 FREQUENCY (HZ) 10* (c) (d) Fig. 4.6: Weighting function A^(to). Zero sequence mode. a) Magnitude, b) Phase angle, c) Real part d) Imaginary part. alternative, however, to having to perform numerical integrations in order to evaluate the function a^(t) and the convolution integrals. This alternative is based on a property of the convolution of an arbitrary function with an exponential function, as discussed next. 61 ID" IO- II- 10" 10" 10" ID" 10* fREOUENCT IHZI 10* .... i MIIII nun 10 I Or iir 10" FREOUENCT IHZI 10* l[f 10" (a) (b) I III II I I llll I I I III II—I I III —I ml II—l l Mill • , ., ., 10" I 10 ID" IO" IO' IO" ID" 10' 10" FREOUENCT IHZI , . I MII II I Mini—I 11II I I I III ll Iff Iff 10 )tt l(f 10 Iff FREOUENCT IHZI (c) (d) F i g . 4.7: W e i g h t i n g f u n c t i o n A ^ ( i o ) . P o s i t i v e sequence mode, a) Magnitude, b) Phase a n g l e , c) R e a l p a r t , d) Imaginary p a r t . 4.3.2 R e c u r s i v e C o n v o l u t i o n I f a t each time s t e p o f a d i s c r e t e p r o c e s s o f s o l u t i o n t h e f o l l o w i n g c o n v o l u t i o n i n t e g r a l has to be e v a l u a t e d : s ( t ) = /°° f ( t - u) K e _ a ( u " T ) d u , T (4.26) 62 then s ( t ) can be d i r e c t l y o b t a i n e d from the known v a l u e o f t h i s f u n c t i o n a t the p r e v i o u s time s t e p , s ( t - A t ) , and the known h i s t o r y o f the f u n c t i o n f at T and (T + At) u n i t s of time e a r l i e r , as f o l l o w s : s ( t ) = ms(t - At) + p f ( t - T) + q f ( t - T - At) . (4.27) In t h i s r e l a t i o n , m, p, and q are c o n s t a n t s t h a t depend on k, a , the i n t e g r a -t i o n s t e p A t , and the n u m e r i c a l i n t e r p o l a t i o n t e c h n i q u e used t o approximate the f u n c t i o n f i n an i n t e r v a l between ( t - T - At) and ( t - T ) . ( T h i s c o n v o l u t i o n p r o p e r t y i s d e r i v e d i n App. I I . 1 . ) . The e v a l u a t i o n o f s ( t ) from r e l a t i o n 4.27 i n s t e a d o f e v a l u a t i n g the complete i n t e g r a l 4.26 (at each time s t e p of the s o l u t i o n p r o c e s s ) can save a c o n s i d e r a b l e amount of computer time. The r e c u r s i v e c o n v o l u t i o n p r o p e r t y i n the m o d e l l i n g o f f r e q u e n c y dependence i s used by Umoto and Hara [ 3 9 ] , by Meyer and Dommel [13] ( t o e v a l u a t e the t a i l p o r t i o n s o f the c o n v o l u t i o n s i n e q u a t i o n s 3.25 and 3.26), and has been s t r o n g l y advocated by Semlyen i n h i s c o n t r i b u t i o n s t o the s o l u t i o n o f t h i s problem (e.g. [ 1 4 ] , [16], and [ 3 0 ] ) . In o r d e r t o a p p l y the p r o p e r t y o f r e c u r s i v e c o n v o l u t i o n t o e v a l u a t e the h i s t o r y i n t e g r a l i n e q u a t i o n 4.22, the f u n c t i o n a ^ ( t ) has f i r s t t o be approximated as a sum of e x p o n e n t i a l terms. However, a sum of e x p o n e n t i a l s i n the time domain c o r r e s p o n d s to a r a t i o n a l f u n c t i o n i n the f r e q u e n c y domain and, t h e r e f o r e , the a p p r o x i m a t i o n can be performed d i r e c t l y on the known f u n c t i o n A^(co) , thus s p a r i n g the need f o r a n u m e r i c a l i n v e r s e t r a n s -f o r m a t i o n t o o b t a i n a ^ ( t ) . B e s i d e s the advantage o f a v o i d i n g t h i s n u m e r i c a l i n v e r s e t r a n s f o r m a t i o n , and due t o the forms o f A^(to) and a ^ ( t ) , the approxima-t i o n o f A^(oi) by a r a t i o n a l f u n c t i o n i s much e a s i e r t o p e r f o r m w i t h a c c u r a c y than the a p p r o x i m a t i o n o f a1(t) by a sum of e x p o n e n t i a l terms. 63 Rational approximation vs. numerical integration: In addition to the savings in computer time, the technique of rational approximation of A^(cu) and recursive evaluation of the convolution integrals can also yield more accurate results than the direct numer-i c a l evaluation of these integrals. In this respect two aspects must be considered. F i r s t , in a direct evaluation of the convolution integrals, the function a^(t) has to be obtained by numerical inverse transformation of A ^ ( i o ) . But, even though in the overall the real and imaginary parts of A^(UJ) in the present formulation are less oscillatory than the corresponding functions in Budner's or Snelson's formulations, the upper part of the frequency range of these functions is s t i l l very oscillatory (e.g. figs. 4.6 and 4.7). This oscillatory behaviour becomes more pronounced when the damping becomes smaller. This i s the case of aerial modes (f i g . 4.7 as compared to f i g . 4.6) and of short lines. Accurate numerical evaluation of the inverse Fourier transform in the very oscillatory regions can be f a i r l y d i f f i c u l t , requiring a large number of frequency points. On the other hand, the magnitude of A - ^ ( c o ) (e.g. figs. 4.6(a) and 4.7(a)) i s very smooth over the entire frequency range, regardless of the line mode or length. The nimerical technique (described in detail later) employed in this work to find a rational approximation for the function A ^ ( c o ) works directly with the magnitude function, whereas the phase function becomes automatically defined from the magnitude function. Since the magnitude function i s evenly smooth, the entire frequency range of A ^ ( c o ) can be evenly approximated and a very accurate closed-form series of exponentials can be obtained for a 1 ( t ) . 6 4 The second aspect of the procedure to determine the history functions (equation 4 . 2 2 ) is the numerical evaluation of the convolution integrals at each time step of the network solution. In this respect, even assuming that a^(t) had been accurately determined from numerical inverse transforma-tion of A^(O J ) , a^(t) would only be available at discrete time points. Because of the sharp spike form of this function ( f i g . 3 . 8 ( a ) ) , the fact of having i t s values at only discrete points may introduce significant additional errors in the evaluation of the convolution integrals. As before, this consideration becomes more important in the case of aerial modes and short lines, where the spike form of a^(t) becomes very high and narrow, tending, in the limit case of no losses, to an impulse. The technique of rational approximation of A-^ (ca) permits a^(t) to be directly obtained in a closed-form;..as a consequence, no additional interpolation errors are introduced in the evaluation of the history functions. Before leaving this point, i t is interesting to mention a possible improvement to the method of direct evaluation of a^(t) by numerical inverse transformation of A^(OJ) . From f i g . 3.8(a), a^(t) can be expressed as a x(t) = p(t - T ) , ( 4 . 2 8 ) where p(t) has the same form as a^(t), but is displaced T units of time towards the origin. From the time shifting property of the Fourier trans-form, the frequency domain form of relation 4 . 2 8 is given by A 1(OJ) = P(oj)e" j a ) T . ( 4 . 2 9 ) The function P(to) is then given by P(u)) = A 1(aj)e j t' 3 T . ( 4 . 3 0 ) 65 The real and imaginary parts of P(co) are much smoother than those of A^(to) over the entire frequency range. Numerical inverse transformation of P(to) would then be much easier to perform with accuracy than that of A^(co). Once p(t) had been determined, the function a^(t) would simply be given by relation 4.28. Further research would be necessary in order to assess the possibilities of this method. The method of rational approximation of A^(co) and recursive eval-uation of the convolution integrals is the one that has been adopted in this work. In the following section the form of the approximating function is considered. The corresponding numerical techniques for obtaining the approximation are discussed in the subsequent section. 4.3.3 Rational Approximation of A-L(OJ) As introduced at the end of the preceding section, the function a^(t) can be expressed as a x(t) = p(t - T ) , (4.31) from where, A 1( U) = P(co)e"jt°T . (4.32) In order to express a^(t) as a sum of exponentials, i t is f i r s t assumed that the function P(s), corresponding to P(co) in the complex plane, can be approximated by a rational function of the form (s + z ) (s + z„) ... (s + z ) P (s) = H . , -\, . , , . " . (4.33) a v (s + P 1)(s + p 2) ... (s + p m) (The subscript "a" i s used to indicate "approximating function.") 66 S i n c e , as shown by e q u a t i o n 3.63 and the c i r c u i t i n f i g . 3.9, the f u n c t i o n A ^ ( t o ) c o r r e s p o n d s to the response of a p a s s i v e , p h y s i c a l system, the f u n c t i o n P (s) i n e q u a t i o n 4.33 must be such t h a t n < m, and the r e a l p a r t a. of the p o l e s must l i e i n the l e f t - h a n d s i d e of the complex p l a n e . Because of n u m e r i c a l c o n s i d e r a t i o n s ( t o be e x p l a i n e d l a t e r ) , the p o l e s and z e r o s are f u r t h e r r e s t r i c t e d to be n e g a t i v e , r e a l , and s i m p l e . Expanding P a ( s ) i n t o p a r t i a l f r a c t i o n s : k, k- k P (s) = + / + ... + . (4.34) a s + p^ s + p 2 s + p m The c o r r e s p o n d i n g time-domain form o f t h i s e q u a t i o n , assuming t h a t the cau-s a l i t y c o n d i t i o n i s s a t i s f i e d (p ( t ) = 0 f o r t < 0 ) , i s g i v e n by a -p-it - p ? t -P t p ( t ) = [k.e + k „ e + . . . + k e m ]U(t) , (4.35) a x z m where U(t) i s the u n i t s t e p f u n c t i o n . From e q u a t i o n s 4.31 and 4.35, and w i t h a ^ a ( t ) a s the f u n c t i o n a p p r o x i m a t i n g a ^ ( t ) , i t then f o l l o w s t h a t -p ( t - x ) -p ( t - x ) -p ( t - x) a, ( t ) = [k,e + k „ e + ... + k e ]U(t - x ) . l a 1 2 m (4.36) T h i s f u n c t i o n has the form r e q u i r e d by e q u a t i o n 4.26 to a l l o w r e c u r s i v e e v a l u a t i o n o f the h i s t o r y c o n v o l u t i o n i n t e g r a l ( e q u a t i o n 4.22). 4.4 N u m e r i c a l S y n t h e s i s of R a t i o n a l A p p r o x i m a t i o n s f o r Z c ( i o ) and A^(OJ ) 4.4.1 G e n e r a l C o n s i d e r a t i o n s A r a t i o n a l f u n c t i o n i n the complex p l a n e has the g e n e r a l f o rm F ( s ) , (4.37) 67 where P(s) and Q(s) are polynomials in the complex variable s = a + joj. If F(s) corresponds to the response of a physical system, for instance, an immittance or transfer function, some restrictions are imposed on the rela-tion between the degrees of P(s) and Q(s) and on the location of the roots (zero and poles) of these polynomials. The problem of approximating a given function of frequency by the general rational function 4.37 is that of finding the polynomials P(s) and Q(s). This problem is not straightforward because, even with the restrictions imposed by the physical system, there i s s t i l l an i n f i n i t e number of possible solutions for P(s) and Q(s). As a result, the range of possible methods of solution is also very large. In general, as in other approximation problems, the question as to which method of solution should be preferred ought to be formulated in terms of obtaining accurate enough results by means of a formula-tion that can be implemented in a practical, "engineering" way with reasonable resources (e.g. costs, man-power, and f a c i l i t i e s ) . Different techniques are mentioned in the classical literature on network synthesis (e.g. Guillemin [28]) for the solution of the rational approximation problem. Examples of these traditional methods are approxima-tions based on the concepts of Butterworth, Chebyshev, Fourier, Taylor, Pad£, and "potential analogies". These methods were developed when limited computa-tional resources were available (before di g i t a l computers) and were mainly applied to the particular problem of the ideal f i l t e r response. The function F(s) in 4.37 was restricted, for example, to P(s) = 1 and to Q(s) being a low order polynomial. Of more recent development (with the development of d i g i t a l computers) are very general optimization techniques [ 3 1 ] , such as, the method of steepest descent, the conjugate gradient method, and the method 68 of Fletcher and Powell. These techniques are based on the optimization of the parameters of a pre-defined F(s) function. Because of the capacity of dig i t a l computers much freedom and generality is allowed to the form of F(s). However, this has created a new problem: too much freedom begets anarchy, and these methods often present convergence problems. The user actually has to supply, somehow, a starting form which is "close enough" to the solution so that the optimization process converges to a "better" fi n a l solution. In many approximation problems, this task is not easy and requires experience and familiarity with both the form of the function being approximated and the internal characteristics of the numerical algorithms. Quoting Szentirmai [32]: "Except for simple cases, the author has been able to coax reasonable results out of iterative optimization programs only by trial-and-error methods of parameter selection." Though this skeptical comment was made by Szentirmai in 1971, i t s t i l l reflects the basic limitation of optimization methods: in many problems, a good guess of the fi n a l form has to be supplied beforehand. In the modelling of frequency dependence in transmission lines, the form of the functions to be approximated depends on the particular line, i t s length, and the particular mode. A pre-selected form of approximating function, or one developed for a particular case, w i l l not in general re-present the best form for other cases, and unless adjustments are made by the user, optimization algorithms may give very poor results or not converge at a l l . Before presenting the technique employed in this work for the rational approximation of the Zc(co) and A^(OJ) functions, reference w i l l be made to approximation techniques employed by other authors in connection with frequency dependence formulations. 69 4.4.2 Approximation Techniques by Other Authors Semlyen makes extensive use of the concept of recursive evaluation of the convolution integrals in his formulations of the frequency-dependence problem. These formulation are in terms of incident and reflected waves and involve the "system propagation function" and the "characteristic admittance function". The approximations used by Semlyen to model these functions have been restricted, however, to a low number of exponential terms (low order Q(s)) in equation 4.37), thus limiting the accuracy to specific frequency regions. In reference [14], the exponentials are obtained from response functions in the time domain; the propagation function is approximated by two exponentials and the characteristic admittance by a constant plus one exponential term (an additional exponential for the characteristic admittance is mentioned in the Closure to reference [14]). Later, in references [16] and [17], the propagation function is approximated directly in the frequency domain by means of a rational function with a Q(s) of order three; the characteristic admittance, however, is assumed to be constant [17]. In Semlyen's most recent publication [30], the need for using higher order approximations for the propagation function (up to six exponentials) is re-cognized as desirable. The method mentioned in this reference to obtain the rational approximations is a modified Gram-Schmidt algorithm (orthogonal polynomials) for a least-square optimization. Problems of ill-conditioning and accuracy in determining the most adequate form of the approximation are also reported in this reference. In another contribution to the modelling of frequency dependence, Ametani [15] suggests a piece-wise straight-line approximation of the system propagation and characteristic impedance functions in the time domain. D i f f i -culties in the use of this method, however, have been indicated by the Bonneville Power Administration; one of the major d i f f i c u l t i e s has been the 70 transition between i n i t i a l sinusoidal steady states and transient states. 4.4.3 Rational Approximation Technique in the New Formulation The method implemented in this work for the synthesis of the equivalent characteristic impedance and for the rational approximation of the history weighting function i s a very flexible and straightforward technique. It is not limited to low-order approximations, as the classical techniques, and does not require the user to supply the adequate parameters and functions, as in modern general optimization techniques. The order and parameters of the approximation are automatically determined by the routine as the approximating function "freely" adapts i t s e l f to the form of the approximated function. Since the form of the approximant is not determined beforehand, numerical instability and accuracy problems due to wrong choices are avoided. Also the entire frequency range of the function can be evenly considered. The technique is based on an adaptation of the very simple concept of asymptotic approximation of the magnitude function, introduced by Bode in 1945 [33], The basic principles for the implementation of this method for the synthesis of Zc(w) and the rational approximation of A^(to) are discussed next. Approximation of the magnitude function: Bode's procedure for approximating the magnitude of a rational function i s illustrated in f i g . 4.8 for the characteristic impedance. The magnitude (in decibels) i s plotted as a function of the logarithm of the frequency. The basic principle of the procedure is to "trace" the approxi-mated curve by straight-line segments. These segments may be either hor-izontal or have a slope that is a multiple of 20 decibels/decade. The 71 points of change of slope (corners or breaking points) define the poles and zeros of the rational approximating function. 5 9 XI 5 7 4 5 6 o a. 0-z 5 4 5 3 5 2 I M i n n j i IN II i i III n — i M i n i — i n u n — i i i n n i i HI n — i i i n u — i MM I I — i M I N I — i i i n n io"5 i r r io"' i io icf itf io" io s io* io 1 io ' FREQUENCY (HZ) Fig. 4.8: Asymptotic approximation of the magnitude of Z (to), zero sec The desired form of Z £^(s) is (from equation 4.12) (s + z )(s + z„) ... (s + z ) z ( s) = H =-— - — eq (s + p )(s + p ) ... (s + p ) (4.38) n Taking the logarithm of the magnitude of this function and multiplying by 20 (to follow the convention of considering the magnitude as gain or attenuation in decibels),it is obtained that 20 log|Z (s)| = 20 logH + 20 log|s + z | + . . . + 20 log|s + z | - 20 log s + p^ - 20 log|s + 20 log s + p (4.39) For s = joj, each one of the terms in this expression has a straight-line asymptotic behaviour with respect to to. Considering, for instance, the term 20 l o g j j u j + Z.J , 72 i t f o l l o w s t h a t l i m f o r co << z^ = 20 l o g z^, t h a t i s , a c o n s t a n t . And, l i m f o r co >> z^ = 20 l o g co, which, as a f u n c t i o n of ( l o g co) , i s a s t r a i g h t l i n e h a v i n g an increment o f 20 d e c i b e l s f o r each u^ /co-^ = 10, t h a t i s , a s l o p e of 20 dB/dec. The combined a s y m p t o t i c b e h a v i o u r of a l l t h e terms i n e q u a t i o n 4.39 can be v i s u a l i z e d as f o l l o w s . Imagining e q u a t i o n 4.38 ( o r e q u a t i o n 4.39) as c o n s t r u c t e d s t e p by s t e p , each time a z e r o c o r n e r ( a t co = z^) i s added, the s l o p e of the asymp-t o t i c curve i s i n c r e a s e d by 20 dB. Each time a p o l e c o r n e r ( a t co = p^) i s added, the s l o p e i s d e c r e a s e d by 20 dB. The f o r e g o i n g d e s c r i p t i o n i s a b r i e f s y n t h e s i s of the well-known fundamentals of Bode's a s y m p t o t i c a p p r o x i m a t i o n . In f i g . 4.8 i t must be n o t e d t h a t the s t r a i g h t - l i n e segments do not r e p r e s e n t the a c t u a l form of the a p p r o x i m a t i n g f u n c t i o n ( e q u a t i o n 4.38) but o n l y an a s y m p t o t i c t r a c i n g of t h i s c u r v e . The a c t u a l a p p r o x i m a t i n g f u n c t i o n i s a smooth c u r v e , w i t h no sharp c o r n e r s , and much c l o s e r to the approximated f u n c t i o n . S i n c e the e n t i r e curve i s t r a c e d d u r i n g the a p p r o x i m a t i o n p r o c e s s , from co -* 0 (the dc c o n d i t i o n can e a s i l y be e x a c t l y matched) to the h i g h e s t f r e q u e n c y at which the approximated f u n c t i o n becomes p r a c t i c a l l y c o n s t a n t , t h e r e a r e no " l e f t o u t " r e g i o n s and a u n i f o r m l y a c c u r a t e a p p r o x i m a t i o n i s o b t a i n e d over t h e e n t i r e f r e q u e n c y range. In the n u m e r i c a l p r o c e d u r e implemented i n t h i s work, the a c c u r a c y o f the f i t t i n g i s f u r t h e r improved by an a l g o r i t h m t h a t s h i f t s the l o c a t i o n of p o l e s and z e r o s about t h e i r f i r s t a s y m p t o t i c p o s i t i o n s . In t h i s p r o c e s s 73 the actual approximating function (equation 4 .38) is compared to the approxima-ted function. The entire process of solution is computationally very fast. The same basic procedure described for the simulation of Z (OJ) is c applied for the simulation of A^(UJ) . Some additional considerations: Since the curves to be approximated (j Z (to) J and |A^(to) | , e.g. figs. 4 .3(a), 4 .4(a), 4 .6(a), and 4.7(a)) are smooth, only real poles and zeros are allowed in the corresponding rational functions (equation 4 .12 for ^ ( to ) and equation 4 .33 for P (to)). Complex poles or zeros can produce local ripples. 3. The problem of local ripples is mentioned by Semlyen [ 3 0 ] , who allows complex exponentials in the approximating functions. Also, in order to obtain an R-C network for the simulation of Z^(to) and simple exponentials in the time domain form of the weighting function,only simple (multiplicity one) poles are permitted. Phase functions: The rational functions 4.12 and 4 . 3 3 , as determined by the method of asymptotic approximation, have no zeros in the right-hand side of the complex plane. Under these conditions, i t is shown in Fourier transform theory (e.g. Papoulis [ 3 4 ] , p. 204) that the corresponding phase functions are uniquely determined from the magnitude functions, and the rational functions belong to the class of minimum phase shift. The agreement between the phases of Z c(to) and P(OJ), and the phases of the corresponding rational approximations obtained in this work confirms the validity of using minimum-phase-shift approximations. The phase angle of the function approximating A^(OJ) is obtained from relation 4 .32 by subtracting -tot from the phase angle of P^(OJ ) : 74 *Ala = *Pa " U T ' ( 4 ' 4 0 ) The time delay T in this relation is directly obtained by comparing the phase angle of P (co) and the phase angle of A (to) : a JL *Pa ~ *A1 ( , . . . x = . (4.41) to If the approximation were exact, that i s , P (to) = P(to) for a l l co, the value 3. of T obtained from relation 4.41 would be the same for a l l frequencies to. In order to minimize the phase error in the approximation, the value of T as given by 4.41 is averaged over the frequency range. Causality condition: The rational approximations Z (s) in equation 4.12 and P (s) in eq a equation 4.33 tend to a constant for s = jco when co -»- 00, and have no poles in the right-hand side of the complex plane. These conditions are enough (e.g. Papoulis [34], pp. 213-214) to assure that the corresponding time domain functions are causal (that i s , the function is equal to zero for t < 0) and thus correspond to the response of a passive physical system. CHAPTER 5 IMPLEMENTATION OF THE SOLUTION 5.1 Recapitulation In the preceding chapters, the synthesis of a new model for the accurate and efficient representation of transmission lines in time-domain transient solutions has been discussed. In the present chapter, the final form of this model and the algorithms for i t s incorporation into Dommel's Electromagnetic Transients Program (EMTP) are considered. In the new model, each line mode (fig . 5.1(a)) is represented by a lumped-parameter equivalent circuit (fig. 5.1(b)) that takes into account the frequency dependence of the line parameters and the distributed nature of the losses. This circuit consists of a constant R-C network that simulates the line characteristic impedance and a voltage source (evaluated at each time step of the system solution) that represents the weighted history of the current and voltage waves travelling along the line. The "weighting" of the travelling signals accounts for the different travelling times and attenuations of the different frequency components. Before proceeding to reduce the line model in f i g . 5.1(b) into i t s f i n a l form, i t is interesting to emphasize the comparison between this circuit and Dommel's simplified one in f i g . 1.3(a). With no losses or frequency depen dence, the characteristic impedance is a simple resistance, that i s , instead of the R-C network in f i g . 5.1(b) the simple resistance R^ in f i g . 1.3(a) is obtained. As for the voltage history source, i t becomes in the simplified cas (equation 3.68) a function of the voltage and current at the other end of the line at only T units of time earlier ("ideal travelling time"); no weighting integrals are thus necessary. 7 6 -*> ' k t vk(t) i m(0<3-l i n e m t (a) N e t w o r k Z e q R 2 (b) Fig. 5.1: (a) Component mode of a transmission line. (b) Equivalent circuit for node k at time step t. After the considerations presented next in this chapter, the line model in f i g . 5.1(b) is further reduced to a single resistance in series with a history voltage source and therefore to exactly the same form as Dommel's simplified •model. The new equivalent resistance "and history source are, however, different from those of Dommel's model and such that the effect of the frequency dependence of the parameters and the distributed nature of the losses is automatically taken into account. 5.2 Evaluation of the Weighted History Source From equation 3.53 , 77 b k ( t ) = { fm ( t ~ u ) a i ( u ) d u • (5.1) But, w i t h a 1 ( t ) approximated as a sum o f e x p o n e n t i a l s (from e q u a t i o n 4.36, changing the n o t a t i o n ) > a l ( t ) = [ k l ' e " B l ( t " T ) + k 2 ' e ^ C t - O + _ + ^ , ^ ( t - T ) ^ _ J } (5.2) I n t e g r a l 5.1 can then be e x p r e s s e d as \ ( t ) = b k ( t ) + b f c ( t ) + ... + b ( t ) , (5.3) 1 2 m where b ( t ) = r f ( t - u) k.' e ' 3 i ( u " T ) d u (5.4) X T f o r i = 1, ..., m. W r i t i n g the c o r r e s p o n d i n g e x p r e s s i o n f o r b ^ a t t h e p r e c e d i n g time i s t e p ( t - At) o f the complete network s o l u t i o n , i t f o l l o w s from the p r o p e r t i e s of r e c u r s i v e c o n v o l u t i o n (App. II.1) t h a t R At- T + At , . b, ( t ) = e ~ P l Z b, ( t - At) + k. / f ( t - u) e P ± < > U " i ; d u . (5.5) X X 1 The o r i g i n a l s e m i - i n f i n i t e i n t e g r a l ( e q u a t i o n 5.4) i s then reduced t o a con-s t a n t times the v a l u e o f the f u n c t i o n a t t h e p r e v i o u s time s t e p p l u s a d e f i n i t e i n t e g r a l o v e r the time s t e p A t . The d e f i n i t e i n t e g r a l i n e q u a t i o n 5.5 can be e v a l u a t e d as f o l l o w s . The f u n c t i o n f i s known from h i s t o r y v a l u e s o n l y a t d i s c r e t e p o i n t s o f time, m u l t i p l e s o f At. S u b s t i t u t i n g f i n e q u a t i o n 5.5 by t h e e q u a t i o n of a s t r a i g h t l i n e p a s s i n g through the known v a l u e s o f f a t [ t - x ] and [ t - (T + A t ) ] and e v a l u a t i n g the r e s u l t i n g c l o s e d - f o r m i n t e g r a l (App. I I . 1 ) , b k ( t ) can be e x p r e s s e d as 78 \ ( t ) = 8 i b k . ( t " A t ) + C i f m ( t ~ T ) + d i f m ( t " 1 " A t> > ( 5- 6) where -3,- At 8 ± = e 1 . (5.7) 1 " 8, and with h. ' = ——;—— , i 3.At ' c i k. ' ~£7 (i - V> <5-8> and k.' d i = - - ^ ( g i - h ! ) . (5.9) Each one of the partial history sources in equation 5.3 can, therefore, be expressed as a function of their corresponding value at the previous time step and the history value of the variables at the other end of the line at ( T) and (x + At) units of time earlier. It is interesting to note that the values of f needed to evaluate the different b, 's are the same: those at m k.£ (t - x) and at (t - x - At). The numerical procedure has included the effect of different travelling times and attenuations in the different coefficients of the partial b^ sources. i 5.3 Equivalent Impedance The equivalent R-C network in f i g . 5.1(b) could be directly i n -corporated in the EMTP, which can directly represent (using the trapezoidal rule of integration) lumped, constant parameters [10]. However, the inter-facing can be simplified by further reducing the R-C network before i t is interconnected with the rest of the system. Two forms of reduction are presented next. The f i r s t one follows directly from the particular nature of the R-C network and leads to an expression that has the same form as the one obtained for the weighted history source. The second reduction is 79 obtained from Dommel's representation of lumped-parameter elements. Direct reduction of the R-C equivalent network: From f i g . 5.1(b), in the frequency domain Ek(co) = I k(u) Z e q(oj) , (5.10) where Z ( c o ) is given by equation 4.12 with s = j c o . The time-domain form eq of equation 5.10 is e k(t) = i k ( t ) * z (t) , (5.11) or, OO e (t) = / i (t - u) z (u)du , (5.12) —oo ^ where z (t) i s the inverse Fourier transform of Z (co) . From equation 4.13 eq eq (changing the notation), and with the causality condition, z (t) = [k n6(t) + k i e " a i t + k 0 e " a 2 t + ... + k e" an t ] i/(t) . (5.13) eq u j. i n With this form for z eq> equation 5.12 can be expressed as e, (t) = e (t) + e (t) + e (t) + . . . + e (t) , (5.14) 0 1 2 n where e k 0 ( t ) - Vk ( t ) - Vk ( t ) ( 5-15 ) and e. (t) = / i . (t - u)k.e a i U d u (5.16) k.. o K. l l for i = 1,2, ..., n. 80 Equation 5.16 has the same form as equation 5.4 with x = 0. The partial voltages in equation 5.14 can thus be expressed as \ 0 " Vk ( t ) (5-17> and e k (t) = m,ev (t - At) + p . i j t ) + q . i j t - At) (5.18) i for i = 1,2, n. As in equations 5.7 to 5.9 the "integration coefficients" are given by m = e i (5.19) 1 - m. and with h. = : , x a.At ' l k. Pi = / ( 1 " V <5-20> 1 and k. q. = - — (m. - h.) . (5.21) x x x Dommel's R-C representation: Using Dommel's [10] representation for resistances and capacitances, each one of the R-C blocks in f i g . 5.1(b) can be modelled as shown in f i g . 5.2. From this model, the relation between e(t) and i(t) can be derived. 81 Fig. 5.2: Dommel's model for R-C combination. This derivation is given in App. II.2. The resulting expression has the same form as equation 5.18, but now the integration coefficients are given by 1 - h (a.At) m i = 1 + H ( c^At) <5'22> and At/2 P i = q i = k i (1 + h (a.Atj ' ( 5' 2 3) Actually, by analogy of the corresponding expressions, the same form of integration coefficients can also be used for the computation of the partial history sources b in equation 5.6. That i s , * i 1 - h (3±At) g i = 1 + H (3^t) ( 5 - 2 4 ) I. I and c. = d. = k.' d ^ .. \) . (5.25) l l l v l + \ (g ^ t ) v ' Both sets of coefficients, the ones obtained from direct integration and the ones derived from Dommel's R-C combination, have been tested in this work. Even though in principle, the direct integration coefficients would be expected to give more accurate results, no significant differences were found 82 i n the t e s t cases s t u d i e d . 5.4 Some G e n e r a l Comments on N u m e r i c a l I n t e g r a t i o n The r e c u r s i v e i n t e g r a t i o n c o e f f i c i e n t s f o r the e v a l u a t i o n of the h i s t o r y s o u r c e b ^ ( t ) ( e q u a t i o n s 5.7 t o 5.9) and f o r t h e e v a l u a t i o n of the v o l t a g e a c r o s s ( e q u a t i o n s 5.19 to 5.21) were o b t a i n e d from l i n e a r i n t e r -p o l a t i o n between the known v a l u e s ( a t d i s c r e t e time s t e p s ) o f the v a r i a b l e f u n c t i o n s f (t) and i ( t ) , and e x a c t e v a l u a t i o n of the c o r r e s p o n d i n g i n t e -g r a t i o n s w i t h the e x p o n e n t i a l f u n c t i o n s (e.g. e q u a t i o n 5.5). Semlyen [30] s u g g e s t s t h a t a second o r d e r i n t e r p o l a t i o n o f the v a r i a b l e f u n c t i o n s s h o u l d l e a d to more a c c u r a t e r e s u l t s . However, i t seems to be a paradox of n u m e r i c a l computations (COSERS [ 3 5 ] , p. 54) t h a t i n many problems the more " r e f i n e d " the d i f f e r e n t i a l s a p p r o x i m a t i n g the d e r i v a t i v e s a r e , the more l i k e l y i t i s to encounter nu-m e r i c a l i n s t a b i l i t y problems. T h i s - sometimes d i f f i c u l t t o e x p l a i n -phenomenon has been encountered many times i n the a n a l y s i s of power systems. As mentioned by Hornbeck ( [ 3 6 ] , p. 165): "One p a r t i c u l a r l y d i f f i c u l t c l a s s o f f u n c t i o n s to i n t e g r a t e n u m e r i c a l l y i s those which are [ s i c ] h i g h l y o s c i l l a t o r y . . . i t s h o u l d be noted t h a t the t r a p e z o i d a l r u l e has c e r t a i n h i g h l y d e s i r a b l e p r o p e r t i e s f o r p e r i o d i c f u n c t i o n s . . . " . These o b s e r v a t i o n s have been c o r r o b o -r a t e d i n p r a c t i c e by the s u c c e s s of the s i m p l e t r a p e z o i d a l r u l e i n p r o v i n g n u m e r i c a l l y s t a b l e and a c c u r a t e f o r t r a n s i e n t s t u d i e s i n power systems (e.g. Dommel, EMTP; Dommel and Sato, t r a n s i e n t s t a b i l i t y [ 3 7 ] ) . The advan-tage of s i m p l i c i t y i n n u m e r i c a l a p p r o x i m a t i o n s has a l s o been found i n the p r e s e n t work i n c o n n e c t i o n w i t h the s i m u l a t i o n of the c h a r a c t e r i s t i c im-pedance and w e i g h t i n g f u n c t i o n through a s y m p t o t i c a p p r o x i m a t i o n s . 83 For the point under consideration at this moment, that i s , the coefficients in the recursive evaluation of the convolution integrals, reported experience in the use of second order interpolation coefficients seems to indicate that numerical instability problems have been encountered in some types of studies. On the other hand, in a l l the cases studied in this work, using linear interpolation coefficients or Dommel1s R-C coefficients, no problems of numerical instability have occurred and the results obtained have been very accurate. In connection with the choice between these last two sets of coefficients, further work is needed in order to establish more de-f i n i t e l y the advantages or disadvantages of one form as compared to the other. 5.5 Reduced Equivalent Circuit With b ^ C t ) expressed by relations 5.3 with 5.6, and ^ ( t ) expressed by relations 5.14 with 5.17 and 5.18, the line equivalent circuit in f i g . 5.1(b) can be simplified as follows. For the voltage across Z : ° eq e k ( t ) = e k 0 ( t > + 1 v t } = R Q i k ( t ) + (Z P ±) i k ( t ) + (Zq±) i k ( t - At) + {Em.e (t - At)} , (5.26) where the summations are from i = 1 to i = n . Combining coefficients and separating present from history quantities, \ ( t ) = V k ( t ) + + W t } > ( 5 - 2 7 ) where 84 and \ = R 0 + ZP± + P = equivalent constant resistance, - k h c ( t ) = ( Z q i ) : L k ( t " A t ) = q \ ( t ' A t ) > = history of current in Z eq •e-khv(t) = Zm-- e^ ( t " A t> (5.28) (5.29) i k. (5.30) = history of partial voltages in Z eq (Underscoring is used to emphasize that the quantities are completely de-fined from past values, that i s , their value at time step t is known.) For the weighted history source b (t ) : K b k ( 0 = Zb k.(t) = f m ( t - x)(Ic i) + f m ( t - x - At) (Ed..) + {Ig.b k.(t - At)} , (5.31) where the summations are from i = 1 to i = m. Combining coefficients, b, (t) = cf (t - x) + df (t - x - At) + Eg.b k.(t - At), (5.32) —k m m I ^ i where f (t - x) = 2v (t - T) - b (t - T) m m m (5.33) and similarly for f (t - x - At) From relations 5.27 and 5.32 the circuit in f i g . 5.1(b) yields the circuit in f i g . 5.3(a) which can be reduced to the simple shunt form shown in f i g . 5.3(b). This reduced shunt model has exactly the same form as ' * ( t j R k §khC<t>*ekhv(t) k O W A ffiV vk(t) ek(t) i k ( t ) vj(t) R k Q ) i k h ( t ) (a) (b) Fig. 5.3: New model reduced equivalent circuit, (a) Series form, (b) Shunt form. 85 Dommel's r e p r e s e n t a t i o n f o r t h e s i m p l i f i e d l i n e ( f i g . 1 . 3 ( b ) ) , and t h e r e f o r e i t can d i r e c t l y be s u b s t i t u t e d f o r t h e s i m p l i f i e d model i n 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. The parameters o f the new l i n e model a r e : = e q u i v a l e n t c o n s t a n t r e s i s t a n c e ( e q u a t i o n 5.28). and ^ ^ ( t ) = t o t a l h i s t o r y s o u r c e - r%hc(t) + % h v ( t ) + V t ) ] / \ (5-34> (e. ( t ) from e q u a t i o n 5.28, e. , ( t ) from e q u a t i o n 5.30, and b. ( t ) from —khc —khv --k e q u a t i o n 5.32). 5.6 A l g o r i t h m of S o l u t i o n The p r o c e s s i n g of frequency-dependent l i n e s w i t h i n the t i m e - s t e p l o o p of the g e n e r a l network s o l u t i o n can be summarized as f o l l o w s : a) B e f o r e s o l v i n g the system network, e v a l u a t e from the r e c o r d e d v a l u e s a t the p r e v i o u s time s t e p s : e k h c ( t ) = q i k ( t - A t ) , (5.35) % h v ( t ) = f m i e k ( t " A t ) ' ( 5 , 3 6 ) i b ^ t ) = ? g.b k ( t - A t ) , (5.37) and 1 b^Ct) = c f m ( t - x) + d f m ( t - T - At) + ^ ( t ) , (5.38) from which, the h i s t o r y c u r r e n t s o u r c e i s g i v e n by i k h ^ " + ^ k h v ( t ) + ^ ( t ) ] / R k • ( 5 - 3 9 ) b) S o l v e the system w i t h the l i n e b e i n g r e p r e s e n t e d by the c i r c u i t i n f i g . 5 . 3 ( b ) . 8 6 c) After solving the system, update the history values needed for the sub-sequent time steps: e k ± ( t ) = (t - At) + P j i j ^ t ) + q 4 i k ( t - At) (5.40) b (t) = g b (t - At) + c.f ( t - T) + d.f (t - T - At) (5.41) fm™ = 2 % ^ ) - V t ) (5.42) The same process is performed for node m. 5.7 I n i t i a l Conditions In the f i r s t steps of the transient study, evaluation of the history current source ( f i g . 5.3(b)) requires the value of the system variables to be known at times before the process is started, that i s , the i n i t i a l conditions at t < 0. The EMTP has a supporting routine to solve the system for i n i t i a l sinusoidal steady-state conditions. Non-sinusoidal conditions have to be supplied by the user. I n i t i a l sinusoidal steady-state conditions for the new line models are evaluated as follows. From the equivalent circuit in the frequency domain in f i g . 4.1(a), the corresponding phasor circuit for sinusoidal steady-state conditions, shown in f i g . 5.4, is directly obtained (upper-scored quantities are phasors). In this circuit V, and I, are known from the i n i t i a l I k Z e q O 1 Fig. 5.4: Equivalent circuit for sinusoidal steady-state conditions at node k. 87 s t e a d y - s t a t e s o l u t i o n of the system g i v e n by the EMTP. F o r the e v a l u a t i o n o f the h i s t o r y c u r r e n t s o u r c e ( e q u a t i o n 5.34) i n the t r a n s i e n t c i r c u i t ( f i g . 5 . 3 ( b ) ) , t h e f o l l o w i n g a d d i t i o n a l q u a n t i t i e s a r e needed ( i n the time and phasor domains): e k . ( t ) - — - V x x b k ( t ) < • B k x i f ( t ) < • F m m In the c i r c u i t i n f i g . 5.4 but from e q u a t i o n 4.13, where Z e q = Z 0 + Z 1 + Z 2 + + Z n ' ( 5 ' 4 4 ) Z0 " k o " R o < 5 - 4 5 > and k "1 jo) + p z i = I „ • (5.46) i E q u a t i o n 5.43 can thus be w r i t t e n as E k " ( Z0 + Z l + Z2 + + V \ " Z 0 I k + Z l I k + - - + Z n V ( 5 - 4 7 ) and from the d e f i n i t i o n of the p a r t i a l v o l t a g e s i n Z ( e q u a t i o n 5.14), V = Z i \ ( 5 ' 4 8 ) x 88 f o r i = 1, . .. , n. The p a r t i a l h i s t o r y s o u r c e s B^ a r e o b t a i n e d as f o l l o w s . 1 From e q u a t i o n 3.48 \ " h V ( 5 ' 4 9> but from e q u a t i o n s 4.32 and 4.34, can be e x p r e s s e d as h = X 1 ( 1 ) + I 1 ( 2 ) + + X l ( m ) ' < 5' 5 0> where k l ( i ) j t o + p ^ — . (5.51) x E q u a t i o n 5.49 can then be w r i t t e n as \ = ( X 1 ( 1 ) + X 1 ( 2 ) + + X l ( m ) ^ m = A.,,. F + A w o . F + ... 4- A, , ,F , (5.52) 1(1) m 1(2) m l(m) m' v J and from the d e f i n i t i o n o f the p a r t i a l h i s t o r y s o u r c e s , e q u a t i o n 5.3, Y - A l ( i ) F m <5'53> x f o r i = 1, ..., m. F can be found from r e l a t i o n 3.42: m F = V + Z I . (5.54) m m eq m From t h e phasor v a l u e s o b t a i n e d i n e q u a t i o n s 5.48, 5.53, and 5.54, the time domain v a l u e s f o r s t e a d y - s t a t e c o n d i t i o n s a t t < 0 a r e g i v e n by e k ( t ) = \\ I c o s < u t + > ' (5.55) i i k. b k ( t ) = l \ I c o s ( w t + > > (5.56) i i k. 89 and f ( t ) = F cos (cut + *_ ) . TTl TTl N ' P r n ' m 1 m (5.57) Analogous r e l a t i o n s a r e o b t a i n e d f o r the c o r r e s p o n d i n g e q u i v a l e n t c i r c u i t at node m. 5.8 I n t e r f a c e w i t h the EMTP The models d e s c r i b e d i n t h i s c h a p t e r have been i n c o r p o r a t e d i n t o the UBC v e r s i o n o f Dommel's e l e c t r o m a g n e t i c t r a n s i e n t s program (EMTP) [38]. T h i s i m p l e m e n t a t i o n i n c l u d e s t h r e e a s p e c t s : p r o c e s s i n g of frequency-dependent l i n e s . G e n e r a t i o n of the l i n e and model parameters : The s t e p s r e q u i r e d t o o b t a i n the l i n e and model parameters a r e as f o l l o w s : a) The parameters R, L, and C f o r the d e c o u p l e d l i n e modes as a f u n c t i o n o f f r e q u e n c y a r e o b t a i n e d from the UBC v e r s i o n of the L i n e C o n s t a n t s Program [ 2 6 ] . The i n p u t d a t a f o r t h i s program a r e t h e c h a r a c t e r i s t i c s o f the c o n d u c t o r s and the tower c o n f i g u r a t i o n o f the a) G e n e r a t i o n o f the l i n e parameters. b) G e n e r a t i o n o f the model parameters. c) M o d i f i c a t i o n of the e x i s t i n g EMPT code f o r the l i n e . b) The d a t a o b t a i n e d i n (a) a r e p r o c e s s e d by FDLINE to o b t a i n 90 c) From the data obtained i n (b) , the c h a r a c t e r i s t i c impedance function Zc(to) i s simulated by TRAN.ZC. The parameters obtained from t h i s simulation: R Q , k^ ., and (see equations 5.17 to 5.21) are stored i n the data f i l e DA.TRAN. d) From the data obtained i n (b) , the weighting function A^(to) i s simulated by TRAN.A1. The parameters obtained from t h i s simulation, k^' and 3^ (see equations 516 to 519) are stored i n the data f i l e DA.TRAN. The parameters i n the data f i l e DA.TRAN are read i n and processed by the new EMTP when the l i n e mode i s i d e n t i f i e d as frequency dependent. Processing of frequency-dependent l i n e s by the EMTP: The e x i s t i n g form of data input for the EMTP has been maintained. A l i n e mode to be simulated according to the frequency-dependence models i s i d e n t i f i e d by spe c i f y i n g MODEL = - 1 i n the corresponding branch card ( f i g . 5.5). The rest of the parameters i n this card correspond, as i n the Model 53 54 F i g . 5.5: Frequency-dependent l i n e mode. e x i s t i n g v ersion, to 60 Hz conditions. For frequency-dependent l i n e s , these parameters are only used f o r i n i t i a l steady-state conditions. When MODEL = - 1 i s detected i n a branch card, the parameters for the frequency-dependence model are read i n from the data f i l e DA.TRAN. 91 The main additions to the existing EMTP code for the simulation of fre-quency-dependent lines using the models developed in this work are the following: a) Evaluation of the coefficients and equivalent resistance. b) Initialization of the history vectors from i n i t i a l conditions. c) Determination of the equivalent current sources. d) Incorporation of the new model parameters into the system network. e) Updating of the history vectors. Some computational aspects: The computer f a c i l i t i e s used in the development of the programs are those of the University of British Columbia. The central unit in the U.B.C. Computing System is an Amdahl 470 V/6 Model II that is run under the MTS operating system. The programs were written in Fortran IV. An approximation of the number of Fortran lines and average running costs required by the frequency dependence routines is given next. The running costs are based on the University rates of $1120 per hour of CPU time under Terminals use, and correspond to the modelling of the zero sequence mode of a typical 100-mi transmission line. For the frequency domain functions 50 points per decade —3 8 are considered in a range from 10 to 10 Hz. a) Frequency-dependent line parameters. Routine name: HWDO:DLINE.0 (UBC Line Constants Program) Running cost: $2.50 b) Computation of Z^Cw) and A^(tu) . Routine name: FLIN Fortran lines: 80 Running cost: $.50 92 c) Simulation of Z (to) . c Routine name: TRAN.ZC Fortran lines: 400 Running cost: $1.00 d) Simulation of A^(OJ) . Routine name: TRAN.A1 Fortran lines: 500 Running cost: $1.50 e) Additions to EMTP code. Fortran lines: 300 Running cost: 10 to 40% more .than with existing code. CHAPTER 6 NUMERICAL RESULTS 6.1 Introduction The results of a series of comparisons and tests for assessing the behaviour of the new line model are presented in this chapter. The results shown in graphical form have been grouped at the end of the chapter and are identified as P.N1.N2 (e.g. P.3.2 refers to plot 2 in Section 3 of this chapter). An index to these graphs is given in pages 122 to 126. 6.2 Frequency-Dependent Line Parameters The transmission line studied in the cases presented in this chapter (with the exception of the field, test comparisons) is a typical 500 kV, 3-phase, single-circuit, transposed, overhead transmission line. The frequency-dependent line parameters are determined by the UBC Line Constants Program. The line input data for this program are: a) Tower configuration: f i g . 6.1 b) Phase conductors: dc resistance = 0.0348 ft/km thickness/diameter = 0.3636 diameter = 40.69 mm c) Ground wires: dc resistance = 1.62 fi/km thickness/diameter = 0.500 diameter = 9.80 mm d) Ground r e s i s t i v i t y = 100 ft-m e) Ground wires are segmented. 94 4t.O> AS. 31' o o o o o o •I o o Fig. 6.1: Tower configuration of reference line. The line constants program uses Carson's formulae [19] to obtain the variation of the phase parameters as a function of frequency and the standard symmetric component transformation to obtain the zero and positive sequence equivalent parameters. The zero and positive sequence series resistance and inductance as a function of frequency were obtained for a range of frequencies between _ o g 10~ Hz and 10 Hz, with 50 frequency points per decade. The variation of these parameters as a function of frequency is shown in graphs P.2.1 and P.2.2. The shunt capacitances and conductances are assumed constant at the following values: zero sequence: C = 7.730 x IO - 3 uF/km G = 0.3 x 10"7 n-1/km positive sequence: C = 14.15 x 10"3 uF/km G = 0.3 x 10 - 7 fi_1/km 6.3 Simulation of the Characteristic Impedance 95 Zero Sequence: The form of the zero sequence characteristic impedance Z c(to) i s shown in the following graphs: i) I Z (u) I -* P. 3.1 i i ) <j>{Z(co)} + P.3.2 i i i ) Re{Z (to) } -»• P.3.3 c iv) Im{Z (to)} P.3.4 c The location of poles and zeros and the reference level of the approximating function Z g q(s) (eqn. 4.12) are shown in table 6.1. The para-meters of the equivalent R-C network (fi g . 4.5) are shown in table 6.2. NUMBER OF ZEROS = 15 ZERO( 1 ) = 0. . 157828E+00 ZERO( 2) = 0 . 308852E + 00 ZERO( 3) = 0. 388782E+00 ZERO( 4) 0. 403171E+00 ZERO( 5) = 0. 404515E+01 ZERO( 6) = 0. 223324E+02 ZERO( 7) = 0. 130467E+03 ZERO( 8) = 0. 541386E+03 ZERO( 9) = 0. 207733E+04 ZERO(10) = 0. 580248E+04 ZERO( 1 1 ) = 0. 124O17E+05 ZERO( 12) = 0. 259181E+05 ZERO(13) = 0. 607560E+05 ZERO(14) = 0. 178462E+06 ZERO( 15) = 0. 119841E+07 NUMBER OF POLES = 15 P0LE( 1) = 0. . 165843E + 00 POLE( 2) = 0. . 3 18982E+00 POLE( 3) = 0. . 4 16563E+00 POLE( 4) = 0. 43G288E+00 POLE( 5) = 0. 375523E+01 POLE( 6) = 0. 209394E+02 POLE( 7) = 0. 122280E+03 POLE( 8) = 0. 507520E+03 POLE( 9) = O. 194768E+04 POLE(10) = 0. 553996E+04 P0LE(11) = O. 120087E+05 P0LE(12) = 0. 250950E+05 P0LE(13) = 0. 588310E+05 P0LEC14) = 0. 172837E+06 P0LE(15) = 0. 115860E+07 MK = 0.455365E+03 Table 6.1: Parameters of Z eq(s). Zero sequence. (Location of corners in Hz, reference level MK in ohms.) 96 ELEMENTS OF FOSTER I EQUIVALENT NETWORK R IN OHMS, C IN MICROFARADS R( 0) 0 455365E+03 C( 0) = ZERO Rl 1) - 276676E+02 C( D - 346858E+05 R( 2) = - 135294E+02 C( 2) = - 368786E+05 R( 3) 0 407864E+02 C( 3) = 0 936750E+04 R( 4) = - 161784E+03 C( 4) = - 225481E+04 R( 5) 0 583433E+02 C( 5) = 0 726428E+03 R( 6) = 0 454623E+02 C( 6) 0 167188E+03 R( 7) 0 431615E+02 c( 7) = 0 301557E+02 R( 8) 0 400909E+02 C( 8) 0 78220GE+01 R( 9) 0 375514E+02 C( 9) = o 217S08E+01 R( 10) = 0 252554E+02 C( 10) = 0 113752E+01 Rl 11) 0 166585E+02 COD 0 795585E+00 R( 12) = 0 1G1149E+02 C( 12) = 0 393555E+00 R( 13) 0 155492E+02 C( 13) = 0 173983E+00 R( 14) = 0 149980E+02 C( 14) 0 613974E-01 R( 15) = 0 155056E+02 C( 15) 0 885923E-02 Table 6.2: Parameters of R-C equivalent network. Zero sequence. The functions Z ( t o ) and Z ( t o ) are compared in the following c eq graphs: i) Magnitudes -> P.3.5 i i ) Phase angles ->- P. 3.6 i i i ) Real parts -* P.3.7 iv) Imaginary parts ->- P.3.8 The magnitude error (in %) is shown in graph P.3.9 and the phase error (in degrees) in graph P.3.10. The maximum magnitude error in most o the curve is 0.2 to 0.3%, with a peak value of about 0.5% in the region o maximum vertical slope. The phase error is also very small, with a peak value of about 0.3° in the region of maximum vertical slope. Positive Sequence: The positive sequence characteristic impedance is shown in the following graphs: 97 i) | z c( u)| i i ) <j>{Z (o>)} c i i i ) Re{Z (to)} c -> P.3.11 -* P.3.12 ->- P.3.13 iv) Im{Z (to)} -> P.3.14 c The location of poles and zeros and the reference level of the approximating function are shown in table 6.3. The parameters R-C of the equivalent network are shown in table 6.4. NUMBER OF ZEROS = 16 NUMBER OF POLES = 16 ZERO( 1) = 0 452926E+00 POLE( 1) = 0 341987E+00 ZERO( 2) = 0 875163E+00 POLE( 2) = 0 813008E+00 ZERO( 3) = 0 128852E+01 POLE( 3) = 0 118584E+01 ZERO( 4) = 0 902754E+00 POLE( 4) = 0 823318E+00 ZERO( 5) = 0 120104E+01 POLE( 5) = 0 109041E+01 ZERO( 6) = 0 156883E+01 POLE( 6) = 0 142411E+01 ZERO( 7) = 0 203965E+01 POLE( 7) = 0 185150E+01 ZERO( 8) = 0 1G7770E+01 POLE( 8) = 0 159506E-"-01 ZERO( 9) = 0 194472E+01 POLE( 9) = 0 185761E+01 ZER0(1O) = 0 225354E+01 POLE( 10) = 0 215220E+01 ZER0(11) = 0 264754E+01 POLE( 11) = 0 253961E+01 ZER0(12) = 0 3153G7E+01 POLE( 12) = 0 302588E+01 ZER0(13) = 0 389754E+01 POLE( 13) = 0 373961E+01 ZERO( 14 ) = 0 42268GE+01 POLE( 14) = 0 499851E+01 ZER0(15) = 0 G54650E+01 POLE( 15) = 0 633909E+01 ZER0(16) = o 242515E+03 POLE( 1G) = 0 239900E+03 MK = 0.237623E+03 Table 6.3: Parameters of Z e q ( s ) . Positive sequence. (Location of corners in Hz, reference level MK in ohms.) ELEMENTS OF FOSTER I EQUIVALENT NETWORK R IN OHMS, C IN MICROFARADS R( 0) = 0 237623E+03 C( 0) = ZERO R( 1) = 0 199350E+03 C( 1) = 0 233450E+04 R( 2) = 0 41 1401E+03 C( 2) = 0 475839E+03 R( 3) = - 589326E+01 C( 3) = - 227740E+05 R( 4) = - 309485E+03 C( 4) = - 624616E+03 R( 5) = O G03876E+02 c( 5) = 0 24 1704E+04 R( 6) = 0 196317E+02 C( G) = 0 569269E+04 R( 7) = 0 142292E+03 c( 7) = 0 604113E+03 R( 8) = 0 241295E+01 C( 8) = 0 413518E+05 R( 9) = - 131995E+03 C( 9) = - 649096E+03 R( 10) = 0 169021E+01 C( 10) = 0 437518E+05 ROD = 0 292379E+01 C(11) = 0 214342E+05 R( 12) • = 0 293638E+01 C( 12) = 0 179125E+05 R(13) • = 0 194440E+01 C( 13) = 0 218881E+05 R(14) = = - 24G582E+02 C( 14) = - 129127E+04 R( 15) = = 0 871910E+01 C( 15) = 0 287953E+04 R( 16) = = 0 258068E+01 C( 16) = 0 257072E+03 Table 6.4: Parameters of R-C equivalent network. Positive sequence. 98 The f u n c t i o n s Z C(OJ) and Z e ^ ( t o ) a r e compared i n the f o l l o w i n g graphs: i ) Magnitudes P. 3.15 i i ) Phase a n g l e s -> P.3.16 i i i ) R e a l p a r t s -* P.3.17 i v ) Imaginary p a r t s -> P.3.18 The o r d e r of the a p p r o x i m a t i o n i s 16. The magnitude e r r o r i s shown i n graph P.3.19 and t h e phase e r r o r i n graph P.3.20. Even though the same t o l e r a n c e , when a s s i g n i n g the asymptotes, was a l l o w e d i n t h i s c ase as i n the case o f the z e r o sequence, the r e s u l t i n g peak magnitude e r r o r i s l a r g e r (about 2.5%). T h i s e r r o r , however, o c c u r s i n the v e r y h i g h s l o p e r e g i o n , w i t h the r e s t o f the cur v e b e i n g more c l o s e l y approximated. A l a r g e r e r r o r i n the h i g h s l o p e r e g i o n i s p r o b a b l y not v e r y s i g n i f i c a n t because o f the f o l l o w i n g c o n s i d e r a t i o n s . Imagining the l i m i t case i n which the s l o p e o f the f u n c t i o n i s a v e r t i c a l l i n e , even a v e r y s m a l l f r e q u e n c y d i s p l a c e m e n t i n the a p p r o x i m a t i n g f u n c t i o n would g i v e an i n f i n i t e e r r o r . In p r a c t i c e , however, the f r e q u e n c y spectrum of the s i g n a l s i s c o n t i n u o u s and a r e l a t i v e l y s m a l l f r e q u e n c y d i s p l a c e m e n t i n the a p p r o x i m a t i n g f u n c t i o n w i l l p r o b a b l y have a v e r y s m a l l e f f e c t i n the f i n a l r e s u l t s . 6.4 A p p r o x i m a t i o n o f the W e i g h t i n g F u n c t i o n The w e i g h t i n g f u n c t i o n A^(co) (eqn. 3.50) depends on the l e n g t h of the l i n e . F o r t h i s r e a s o n , d i f f e r e n t l e n g t h s were c o n s i d e r e d f o r the t e s t i n g o f the f i t t i n g r o u t i n e s . The r e s u l t i n g a p p r o x i m a t i o n s f o r the zero and p o s i t i v e sequence modes of these l i n e l e n g t h s a r e shown i n the f o l l o w i n g graphs: 99 Zero Sequence: i ) 5-mi P .4 .1 i i ) 30-mi P 4 .3 i i i ) 100-mi P 4 5 i v ) 500-mi P 4 7 v) 1000-mi P 4 9 Sequence: i ) 5-mi -»- P 4 2 i i ) 30-mi -»• P 4 4 i i i ) 100-mi -»• P 4 6 i v ) 500-mi P 4 8 v) 1000-mi P 4 10 As shown by these graphs, the a p p r o x i m a t i o n i s v e r y c l o s e t o the f u n c t i o n i n a l l t h e c a s e s . The a p p r o x i m a t i o n s were performed so as to a l l o w a maximum magnitude e r r o r of 0.5% i n the r e g i o n o f the cur v e s between 0 Hz and the fr e q u e n c y a t which the magnitude reaches 0.7. T h i s e r r o r r e g i o n was d e f i n e d so as not t o i n c l u d e the h i g h - s l o p e zone i n the maximum e r r o r c r i t e r i a (even though the r e s u l t i n g e r r o r i n the h i g h - s l o p e zone i s a l s o v e r y s m a l l ) . The number o f p o l e s and z e r o s o f t h e a p p r o x i m a t i n g r a t i o n a l f u n c t i o n s (eqn. 4.33) i s shown i n t a b l e 6.5. l e n g t h 5 mi 30 mi 100 mi 500 mi 1000 mi Seq. z e r o pos. z e r o pos. z e r o pos. z e r o pos. z e r o pos. z e r o s 12 13 11 11 12 13 15 20 17 23 p o l e s 16 14 15 14 20 17 21 25 23 29 T a b l e 6.5: R a t i o n a l a p p r o x i m a t i o n s f o r d i f f e r e n t l i n e l e n g t h s . 100 A more d e t a i l e d study comparing the function and the approximation for a l i n e length of 100 miles i s presented next. Zero Sequence (100 mi) The form of A^(to) for the zero sequence mode and a length of 100 miles i s shown i n the following graphs: i ) |A1(a))| P.4.11 i i ) 4>{A1(co)} P.4.12 i i i ) Re{A (u>)} P. 4.13 iv) Im{A (OJ) } ->• P.4.14 The parameters of the r a t i o n a l approximating function Pa(s) (eqn. 4.33) and the phase displacement T (eqns. 4.31 and 4.32) are shown i n table 6.6. NUMBER OF ZEROS = 15 NUMBER OF POLES = 20 ZER0( 1) = 0 367679E+G2 POLE 1 ) = 0 359264E+02 ZER0( 2) = 0 155040E+03 POLE 2) = 0 151515E+03 ZER0( 3) = O 2G9450E+03 POLE , 3) = 0 263346E+03 ZER0( 4) = O 385017E+03 POLE 4) = 0 376262E+03 ZER0( 5) = O 4959G1E+03 POLE 5) = 0 484659E+03 ZERO( G) = 0 610179E+03 POLE G) = 0 596240E+03 ZERO( 7) = 0 733576E+03 POLE 7) = 0 716917E+03 ZERO( 8) = O 842268E+03 POLE 8) = 0 823059E+03 ZER0( 9) = 0 966982E+03 POLE 9) = o 944968E+03 ZER0(10) = 0 108511E+04 POLE 10) = 0 106042E+04 ZER0(11) = 0 134982E+04 POLE 11) = 0 124524E+04 ZER0(12) = O 184181E+04 POLE( 12) = 0 167985E+04 ZER0(13) = o 239982E+04 POLE 13) = 0 216360E+04 ZER0(14) = O 213423E+04 POLE! 14) = 0 183774E+04 ZER0C15) = O 251200E+04 POLE ( 15) = 0 239400E+04 POLE! 16) * 0 326695E+04 POLE( 17) = 0 119650E+05 P0LE( 18) - 0 712960E+04 = 0.487557E+20 POLE ( 19) = 0 156694E+05 ) = 0.5976E -03 P0LE( 20) = 0 226480E+05 Table 6.6: Parameters of the r a t i o n a l approximation of A^(a)). Zero sequence mode, 100 miles. (Location of corners i n Hz, reference l e v e l MK i n p.u., x i n sec.) In connection with the form of A^(co), the observation made at the end of Section 4.3.2 i n the sense that the o s c i l l a t o r y behaviour of the phase of t h i s function could be circumvented by considering the function P(to) (eqn. 101 4.30) instead of A^(to) i s corroborated by graph P.4.15 which shows the smooth phase variation of P (to) . a The function A^(OJ) and i t s rational approximation are compared in the following graphs: i) Magnitudes -*• P. 4.16 i i ) Phase angles -*• P.4.17 i i i ) Real parts -> P.4.18 iv) Imaginary parts -+ P.4.19 The magnitude error is shown in graph P.4.20 and the phase error in graph P.4.21. As mentioned before, the magnitude error is less than 0.5% in most of the curve. The highest peak error is about 2% in a region of very high vertical slope. For the reasons explained before, however, the higher error in this zone is not very significant, and as can be seen from graph P.4.16,the two curves are practically identical. The phase error along the curve is also very small, of about 0.3°. In the region where the magnitude is practically zero, the magnitude and phase errors lose their significance. Positive Sequence (100 mi): The function A^(to) for the positive sequence mode and a length of 100 miles i s shown in the following graphs: i) |A1(to)| +P.4.22 i i ) cj){A1(w)} •* P.4.23 i i i ) Re{A 1(u)} + P.4.24 iv) Im{A (u>) } -»• P.4.25 The phase angle of the function P ( t o ) ( = A ( t o ) e^ 1) is shown in a J_a graph P.4.26. The parameters of the corresponding rational approximation are shown in table 6.7. 102 NUMBER OF ZEROS = 13 NUMBER OF POLES = 17 ZERO( 1) = 0 923471E+03 P0LE( D = o 902475E+03 ZERO( 2) = 0 852040E+04 POLE( 2) = 0 832608E+04 ZERO( 3) = 0 173949E+05 POLE( 3) = 0 170000E+05 ZERO( 4) = 0 260346E+05 POLE( 4) = 0 254396E+05 ZERO( 5) = 0 347048E+05 POLE( 5) - 0 339150E+05 ZERO( 6) = 0 431969E+05 P0LE( 6) = 0 422 154E+05 ZERO( 7) = 0 519331E+05 POLE( 7) = 0 507520E+05 ZERO( 8) = 0 603169E+05 POLE( 8) = 0 589390E+05 ZER0( 9) = o 684622E+05 P0LE( 9) = o 669031E+05 ZERO( 10) = 0 768163E+05 P0LE( 10) - 0 750654E+05 ZERO(11) = 0 879200E+05 POLE( 11) = 0 783631E+05 ZERO(12) = 0 136188E+0G POLE( 12) = 0 117248E+06 ZERO(13) = 0 312719E+06 POLE( 13) = 0 169485E+06 POLE( 14) = 0 367547E+06 POLE( 15) = 0 512864E+06 POLE( 16) = 0 104240E+07 MK = 0.149967E+24 POLE( 17) = 0 231754E+07 TAU0 = 0 .5394E-03 Table 6.7: Parameters of the rational approximation of A (u>) . Positive sequence mode, 100 miles. (Location of corners in Hz, reference level MK in p.u., T in sec.) The comparisons between the function and its approximations are shown in the following graphs: i) Magnitudes ->- P .4 .27 i i ) Phase angles ->• P .4 .28 i i i ) Real parts ->• P .4 .29 iv) Imaginary parts -»• P .4 .30 The magnitude error is shown in graph P.4.31 and the phase error in graph P.4.32. The apparent increase in the phase error in the upper range of gra P.4.32 is probably due to the lack of enough frequency points to follow the very fast variations of the phase function in this region of the curve. Time domain form of the weighting functions: The time domain form of the zero and positive sequence weighting functions, as obtained from eqn. 4.36 from inverse transformation of the 103 corresponding rational approximations, is shown in graphs P.4.33 and P.4.34. The phase displacement x is 0.5976 msec for the zero sequence mode and 0.5394 msec for the positive sequence mode. Some interesting conclusions can be derived by comparing the zero and positive sequence forms of a^(t). When the functions are compared on the time scale of graphs P.4.33 and P.4.34(a), a^-pos. looks almost like an im-pulse, 6(t-x), whereas the spike form of a^-zero is clearly apparent. This • means that taking, for instance, a time step At for the transient process of about 1/10 of x, say 0.05 msec, the "weighting" effect of a^-zero w i l l be well taken into account, whereas a^-pos. w i l l be almost an impulse. This is equivalent to having A^(to)-pos. constant in the frequency domain, which means that the larger damping of the higher frequency components is not taken into account. Actually, a period of 0.05 msec corresponds to a frequency of 20 kHz, at which point A^(to)-pos. (graph P.4.22) is just beginning to decay while A^ (to)-zero (graph P.4.11) has decayed almost completely. The expanded time-axis plot of a^(t)-pos. in graph P.4.34(b) shows that taking, for instance, a At of about 0.001 msec the weighting effect of this function w i l l be clearly represented. This At corresponds to a frequency of 10^ Hz; at this frequency, A^(to)-pos. (graph P.4.22) has decayed almost completely. 6.5 Analytical Tests: Frequency Domain Short circuit and open circuit terminations constitute extreme loading conditions in a transmission line. These limit conditions provide a simple and yet demanding way of assessing the effect of the approximations for the characteristic impedance and weighting function in simulating the line response in the frequency domain. 1 0 4 From the general solution of the line equations in the frequency domain (eqns. 1 . 3 4 and 1 . 3 5 ) , with a purely sinusoidal (single frequency) voltage source (E g) connected at the sending end of the line,and the receiving end short-circuited (Vm = 0), the current at the sending end in terms of A^OJ) (eqn. 3.50) and z c ( w ) i s given by (in phasorial form): E (to) 1 + A2(co) k Z c ( u ) ) 1 - A2(co) Similarly, with the receiving end of the line open-circuited, the voltage at this end i s given by: 2 A (OJ) V (OJ) = E (OJ) i= . ( 6 . 2 ) 1 + A^(OJ) It is interesting to observe from eqns. 6 . 1 and 6 . 2 that the open-circuit voltage is independent of the characteristic impedance, whereas the short-circuit response i s directly proportional to 1 / Z c (O J ) . This would explain why some frequency dependence models that assume a constant characteristic impedance can give acceptable results when tested only for open-circuit or close-to-open-circuit conditions. On the other hand, the correct modelling of the frequency depedence of Z^(OJ) is very important for short-circuit or near-to-short-circuit conditions. From eqns. 6 . 1 and 6 . 2 , the short-circuit and open-circuit responses of the line using the exact values of Zc(co) and Z^(OJ) can be compared with the responses obtained by using the values of the corresponding approximating functions. These comparisons can be made over the entire frequency range. Also, in order to il l u s t r a t e the importance of taking into account the frequency dependence of the parameters, the short-circuit and open-circuit responses can be calculated assuming constant line parameters. 105 In the results presented next, the following responses were compared: a) Response with exact parameters (C) vs. response with approximated parameters (P). b) Response with exact parameters (C) vs. response with 60 Hz constant parameters (P60). c) Response with exact parameters (C) vs. response with 750 Hz constant parameters (P750). In most transient simulations using constant parameters, these parameters are taken at 60 Hz in order to correctly reproduce the fundamental, 60-Hz signal injected by the generators. In an attempt to more closely simulate the higher harmonics, sometimes the constant parameters are taken at a higher frequency (e.g. 500 to 1000 Hz). As i t w i l l be shown in the comparisons with 750 Hz parameters, however, taking the parameters at these higher frequencies can represent a considerable detriment of the 60-Hz component, without major improvement in the frequency range outside the vicinity of the frequency at which the parameters are taken. Zero Sequence Responses: A) Short-circuit responses: Assuming a source voltage of 100 kV, the magnitudes of the short-circuit currents (in amperes) at node k are compared in the following graphs: i) Exact vs. approximation -> P.5.1 i i ) Exact vs. 60-Hz ->• P.5.2 i i i ) Exact vs. 750-Hz -> P.5.3 As can be seen from these plots, the agreement between the exact and approximated responses is very good. Sixty-hertz constant parameters give 106 a good agreement o n l y from about 10 to 500 Hz, whereas the 750-Hz parameters c o v e r ( l e s s w e l l ) a range from about 100 t o 1000 Hz and g i v e an e r r o r o f about 40% a t the fundamental 60-Hz f r e q u e n c y . The magnitude e r r o r s a r e shown i n the f o l l o w i n g p l o t s : i ) E x a c t v s . a p p r o x i m a t i o n -> P.5.4 i i ) E x a c t v s . 60-Hz -* P.5.5 i i i ) E x a c t v s . 750-Hz P.5.6 As shown i n graphs P.5.4 the maximum e r r o r o b t a i n e d w i t h the a p p r o x i m a t i n g f u n c t i o n s i s l e s s than about 1% i n most of the f r e q u e n c y range. The e r r o r i s r e l a t i v e l y l a r g e r a t about 300 Hz, where, however, i t does not have much s i g n i f i c a n c e s i n c e the magnitude of the c u r r e n t i s r e l a t i v e l y c l o s e to z e r o . Another p o i n t where the e r r o r i s s l i g h t l y l a r g e r than 1% i s a t about 600 Hz which c o r r e s p o n d s to a r e l a t i v e l y narrow peak of the response. The e r r o r s a t v e r y narrow peaks have l e s s s i g n i f i c a n c e than i n the f l a t t e r r e g i o n s s i n c e these p o i n t s a r e almost s i n g u l a r i t i e s of the response f u n c t i o n . The e r r o r i s a l s o r e l a t i v e l y l a r g e r i n the range from about 0.1 to 5 Hz. T h i s e r r o r i n the s h o r t - c i r c u i t response a t r e l a t i v e l y low f r e q u e n c i e s i s v e r y d i f f i c u l t t o c o n t r o l . T h i s i s due to the n u m e r i c a l s e n s i t i v i t y of the -2 -2 f a c t o r (1 + A ^ ) / ( l - A^) i n eqn. 6.1 when A^(to) i s v e r y c l o s e t o 1. F o r i n s t a n c e , f o r A^ e x a c t = 0.9998 and A^ approx. = 0.9999, t h a t i s , a d i f f e r e n c e o f o n l y + 0.01% i n A^, the e r r o r i n the f a c t o r i s + 100%. ( A c t u a l l y , as can be seen i n graph P . 5 . 1 ( a ) , the " e x a c t " response i t s e l f p r e s e n t s a s l i g h t l y e r r a t i c b e h a v i o u r i n t h i s r e g i o n . ) To m i t i g a t e t h i s problem i n the s i m u l a t i o n of the low f r e q u e n c y r e g i o n , the f u n c t i o n A-j_(a))-approx. i s matched w i t h -14 A^(w)-exact a t ID = 0 (dc c o n d i t i o n ) w i t h an e r r o r o f l e s s than 10 %. T h i s r e s u l t s i n a p r a c t i c a l l y e x a c t matching o f the dc s h o r t - c i r c u i t c u r r e n t and 107 helps in controlling the error over the low frequency range. This exact matching of the dc condition i s very important because i t assures the correct simulation of a dc level in the signals. In the case of the constant-parameter approximations (graphs P.5.5(a) and P.5.6(a)) the error at low frequencies is very large and dc levels cannot be simulated. It can also be seen from the error graphs that in the case of the approximation the response at 60 Hz is also exactly matched, thus assuring the correct simulation of the sources frequency. B) Open-circuit responses: Assuming a source voltage of 1.0 (per unit), the open-circuit voltage responses at the receiving end of the line are compared in the following graphs: i) Exact vs. approximation ->- P.5. 7 i i ) Exact vs. 60-Hz -> P.5.8 i i i ) Exact vs. 750-Hz -> P.5.9 As shown in these plots, the approximation response correctly matches the exact response over the entire frequency range. The constant-parameter representations, on the other hand, give a good approximation only for the range between 0 and about 200 Hz. The magnitude errors are shown in the following graphs: i) Exact vs. approximation ->- P.5.10 i i ) Exact vs. 60-Hz ->• P.5.11 i i i ) Exact vs. 750-Hz -»- P.5.12 The error for the approximation i s less than about 1.5% over the entire frequency range, with slight increases in the regions of very high 108 v e r t i c a l slope. The response at the fundamental 60 Hz frequency i s exactly matched. P o s i t i v e Sequence Responses: A) S h o r t - c i r c u i t responses: For a source voltage of 100 kV, the magnitudes of the s h o r t - c i r c u i t currents using the exact, approximating, and constant-parameter Z c(w) and A^(w) functions are shown i n the following graphs: i ) Exact vs. approximation -> P.5.13 i i ) Exact vs. 60-Hz P. 5.14 i i i ) Exact vs. 750-Hz •+ P.5.15 Considering f i r s t the frequency region of graphs (a) (10 t o 10 Hz) i t can be seen that the approximation follows the exact response very c l o s e l y despite the very narrow resonance spikes. The constant-parameter responses, i n contrast with the zero sequence case, also give a r e l a t i v e l y f a i r approxima-4 tion (with the exception of the very narrow spike regions) up to about 10 Hz. Comparing now the low frequency responses (graphs ( b ) ) , i t can be seen that with the exception of the range from about 0.1 to 5 Hz, the approximation follows the exact response very c l o s e l y . As i n the case of the zero sequence, however, the range between 0.1 and 5 Hz presents a larger e r r o r . For the p o s i t i v e sequence the r e l a t i v e l y low frequency range i s even more s e n s i t i v e to very small differences i n A^(co) than f o r the zero sequence because the values of A^(w) are closer to 1.0. It i s i n t e r e s t i n g to observe that for t h i s l i n e length the 60 Hz constant-parameter representation gives a very good simulation of the low-frequency range. This i s not so, however, for the 750 Hz simulation. 109 The comparison between the magnitude errors are shown in the following graphs: i) Exact vs. approximation -> P.5.16 i i ) Exact vs. 60-Hz -*- P.5.17 i i i ) Exact vs. 750-Hz •+ P.5.18 B) Open-circuit responses: Then open-circuit responses for a source voltage of 1.0 p.u. are compared in the following graphs: i) Exact vs. approximation •+ P.5.19 i i ) Exact vs. 60-Hz ->- P.5.20 i i i ) Exact vs. 750-Hz -»- P.5.21 The corresponding magnitude errors are shown in graphs: i) Exact vs. approximation -> P.5.22 i i ) Exact vs. 60-Hz -> P.5.23 i i i ) Exact vs. 750-Hz P.5.24 As i t can be seen from these graphs the approximation response is very good over the entire frequency range, whereas the constant-parameter 4 responses give fair results only up to about 10 Hz. 6.6 Analytical Tests: Time Domain The accuracy of the frequency domain approximations of the characteristic impedance and weighting function in simulating the response of the system was verified in the previous section for short-circuit and open-circuit frequency-response tests. Short-circuit and open-circuit single-frequency conditions also provide a simple way of verifying the numerical 110 behaviour of the time-domain line models. For this purpose, the line, represented by i t s frequency-dependence transient model was energized by a purely sinusoidal voltage source. Starting from sinusoidal steady-state conditions, transient simulations (without disturbances), with the receiving end of the line short- and open-circuited, were run with the EMTP. Since no disturbances were introduced, the time-domain solutions obtained with the EMTP should be (ideally) perfectly sinusoidal waves. The magnitude and phase of these waves should agree with the values obtained analytically from eqns. 6.1 and 6.2 . In order to assess the accuracy of the numerical discretizations involved in the time-domain models, the results of the transient simulations were compared not with the analytical results using exact parameters, but with the analytical results using the values of Z c ( io ) and A^(w) given by the approximating functions. (The differences between the analytical values obtained from the exact and approximate parameters were presented in the previous section.) These time-domain tests were performed for the zero and positive sequence modes of 100 miles of the reference line for different frequencies of the applied voltage source. In a l l the cases studied, the time-domain functions had the correct sinusoidal wave-form,and the magnitude peaks were in very good agreement with the analytical values. Some of the results of these simulations (for the short circuit case) are presented next. In these cases, the simulations were also performed using the EMTP code with no frequency dependence (constant-parameter model). In order to compare the results of the constant-parameter model with the same reference value as in the case of the frequency-dependence model, the line parameters at the source frequency for the constant-parameter model were I l l calculated from the approximated Zc(w) and A^(CJ) functions. Since only one frequency is involved in these simulations, and the parameters at this frequency are used for the constant-parameter model, the results obtained with this model would also be expected to coincide with the analytical values. That this was not exactly the case, however, w i l l be seen from the obtained results. Zero Sequence Comparisons: i) f = 60.2 6 Hz. The peak magnitudes of the simulations are shown in table 6.8. In this table, I ^ = simulation with frequency-dependence model. I = simulation with constant-parameter model, cp I = analytical value from eqn. 6.1. ex J AI,.. and AI = % variation with respect to I fd cp ex At = time step in the transient simulation. t(msec) I f d ( A ) AI f d(%) cp AI (%) cp 3.80 446.8 -0.04 447.3 0.07 12.1 -447.0 0.00 -447.6 0.13 20.4 446.8 -0.04 447.4 0.09 28.6 -447.0 0.00 -447.5 0.11 36.9 446.8 -0.04 447.5 0.11 45.2 -447.0 0.00 -447.5 0.11 Table 6.8: Short ci r c u i t , zero seq., f=60.26 Hz, 1^=447.0 A, At=0.1 msec. As shown in table 6.8 the results of the time domain simulations agree very well with the expected analytical value. i i ) f = 2089 Hz. The results for this case are shown in table 6.9. 112 t(msec) I f d ( A ) A I f d ( % ) cp AI (%) cp .23 -281.0 -0.35 -382.5 36 .47 281.0 -0.35 382.5 36 .71 -281.1 -0.32 -382.5 36 .95 281.1 -0.32 382.5 36 1.19 -281.1 -0.32 -382.5 36 1.42 281.4 -0.21 382.5 36 Table 6.9: Short c i r c u i t , zero seq., f=2089 Hz, I =282.0, At=0.005 msec. ex The large error i n the case of the constant-parameter representation at t h i s frequency i s probably due to the fact that t h i s model does not take into account the d i s t r i b u t e d nature of the system losses. i i i ) f = 1 Hz. The r e s u l t s f o r t h i s case are shown i n table 6.10. t(msec) I f d ( A ) A I f d ( % ) cp AI (%) cp 175 21 058 4.14 21 823 7.93 675 -19 996 -1.11 -21 823 7.93 1175 20 287 0.33 21 823 7.93 Table 6.10: Short c i r c u i t , zero seq., f=l Hz, I =20 220 A, At=0.1 msec. -l ' e x The i n i t i a l o f f s e t i n the case of the frequency-dependence simula-t i o n i s probably due to the l i m i t e d number of d i g i t s (four) that can be entered as input parameters f or the pre-transient steady-state c a l c u l a t i o n s in the EMTP. This i s probably also the reason for the error i n the constant-parameter simulation. No i n i t i a l o f f s e t occurs for t h i s model, however, since the parameters used f o r the steady-state and transient simulations are the same. The l i m i t a t i o n i n the number of d i g i t s f o r the steady-state parameters can be avoided i n the frequency dependence version of the EMTP by c a l c u l a t i n g these 113 parameters internally from the transient parameters. This feature w i l l be implemented in a future version of the program. Positive Sequence Comparisons: i) f = 60.26 Hz. The results of these simulations are shown in table 6.11. t(msec) I f d(A) M f d(%) cp AI (%) cp 4.1 1982 -0.10 1974 -0.50 12.4 -1985 0.05 -1977 -0.35 20.7 1982 -0.10 1974 -0.50 29.0 -1985 0.05 -197 6 -0.40 37.3 1983 -0.05 1974 -0.50 45.5 -1985 0.05 -1976 -0.40 Table 6.11: Short circuit, pos. seq., f=60.26 Hz, I =1984, At=0.1 msec. ex i i ) f = 10 000 Hz. The results for this case are shown in table 6.12. t(msec) I f d(A) AI f d(%) I (A) cp AI (%) cp .027 -771.6 -0.48 -Ilk A -0.12 .077 771.6 -0.48 llk.7, -0.13 .127 -771.6 -0.48 -774,3 -0.13 .177 771.6 -0.48 774.4 -0.12 .227 -771.6 -0.48 -774.5 -0.10 .277 771.6 -0.48 774.6 -0.09 Table 6.12: Short circuit, pos. seq., f=10 000 Hz, 1^=775.3, At=0.001 msec. The constant-parameter, lumped-resistance model gives much better results at high frequencies for the case of the positive sequence mode than i t does for the zero sequence mode. This is probably due to the fact that 114 a t t h e s e f r e q u e n c i e s t h e r e s i s t a n c e as compared to the r e a c t a n c e i s much s m a l l e r f o r t h e p o s i t i v e sequence mode than i t i s f o r t h e z e r o sequence mode. As a consequence, the e f f e c t o f the r e s i s t a n c e b e i n g r e p r e s e n t e d as c o n c e n t r a t e d i n s t e a d o f d i s t r i b u t e d has much l e s s e f f e c t i n t h e s i m u l a t i o n . i i i ) f = 1 Hz. The r e s u l t s f o r t h i s c a s e a r e shown i n t a b l e 6.13. t(msec) I f d ( A ) A I f d ( % ) cp AI (%) cp 70 57 860 -0.50 57 600 -0.95 57 0 -58 130 -0.03 -57 600 -0.95 1070 58 160 0.02 57 600 -0.95 T a b l e 6.13: S h o r t c i r c u i t , pos. seq., f = l Hz, I =58 150 A, At=0.1 msec. 6.7 Examples of T r a n s i e n t S i m u l a t i o n s The importance of c o r r e c t l y s i m u l a t i n g the e n t i r e f r e q u e n c y range of the s i g n a l s (from dc (0 Hz) to t h e v e r y h i g h f r e q u e n c i e s ) i s i l l u s t r a t e d i n t h i s s e c t i o n w i t h some comparisons of t r a n s i e n t s i m u l a t i o n s u s i n g t h e frequency-dependence models and u s i n g c o n s t a n t - p a r a m e t e r r e p r e s e n t a t i o n s at 60 and 750 Hz. In t h e s e examples, the z e r o sequence mode of 100 m i l e s of the r e f e r e n c e l i n e was d i r e c t l y e n e r g i z e d w i t h an i d e a l v o l t a g e source at the s e n d i n g end and v a r i o u s c o n d i t i o n s were c o n s i d e r e d a t t h e r e c e i v i n g end. S h o r t c i r c u i t s i m u l a t i o n s : W i t h z e r o i n i t i a l c o n d i t i o n s , a 60-Hz, 100 kV peak, s i n u s o i d a l v o l t a g e s o u r c e was a p p l i e d t o t h e sending end of the l i n e w i t h t h e r e c e i v i n g end s h o r t - c i r c u i t e d . The a p p l i e d v o l t a g e p a s s e s t h r o u g h z e r o a t t = 0. The r e s u l t i n g s h o r t c i r c u i t c u r r e n t s a t t h e r e c e i v i n g end a r e shown i n the 115 following graphs: i) Frequency-dependence model ->- P.7.1 i i ) Constant 60-Hz-parameter model ->• P.7.2 i i i ) Constant 750-Hz-parameter model P.7.3 Due to the extended time range of these simulations, a relatively large time step, At = 0.4 msec was used. (The travelling time is i = 0.5976 msec.) A most interesting element in these graphs is the simulations of the exponentially decaying dc component of the short-circuit current. This component dies very fast in the 60-Hz model and almost does not appear at a l l in the 750-Hz model. These results corroborate the predicted behaviour of the models from the short circuit frequency responses (graphs P5.1(a), P.5.2(a), and P.5.3(a)). They also show the importance of the line model to be capable of correctly simulating dc conditions. A closer view of the f i r s t cycles of this simulation is shown in the graphs: i) Frequency dependence (P) vs. 60 Hz (P60) •> P.7.4 i i ) Frequency dependence (P) vs. 750 Hz (P750) ->• P.7.5 These simulations were run with At = 0.1 msec. The peak value of the short-circuit current given by the different simulations was, approximately: i) Frequency dependence: 860 A i i ) 60 Hz parameters: 7 65 A (-11%) i i i ) 750 Hz parameters: 280 A (-67%) Another important value is the f i r s t zero crossing (for fast circuit breakers). This value was approximately: 116 1) Frequency dependence: 14.6 msec i i ) 60-Hz parameters: 13.5 msec i i i ) 750-Hz parameters: 9.8 msec A short c i r c u i t simulation for the case when the voltage i s at i t s peak value at t = 0 i s shown in the graphs: i) Frequency dependence vs. 60 Hz P. 7. 6 i i ) Frequency dependence vs. 750 Hz -*• P.7 .7 A At of 0.01 msec was used for these simulations. The differ e n c e i n the harmonic content of the signals for the d i f f e r e n t models i s quite noticeable. Since very l i t t l e dc o f f s e t i s present i n t h i s case, the 60-Hz simulation follows the general shape of the frequency-dependence simulation; the 750-Hz simulation, however, gives very inaccurate r e s u l t s in both magnitude and phase angle. As a t h i r d example of short c i r c u i t modelling, the l i n e was energized with a step voltage of 100 kV. The corresponding receiving end currents are shown i n the following graphs: i) Frequency dependence vs. 60 Hz •*• P.7.8 i i ) Frequency dependence vs. 750 Hz P.7.9 (At-0.01 msec). The step response simulation for the frequency-dependence model was extended over a longer time (graph P.7.10; At = 0.4 msec) i n order to v e r i f y that the model c o r r e c t l y reproduces the dc response. Neglecting the very small shunt current, the f i n a l asymptotic value of the step response should coincide with the value obtained i n the frequency response curves (graph P.5.1(a)) and with the value calculated a n a l y t i c a l l y . Reading the approximated values from the graphs, 117 from graph P.7.10, 1^ = 53 500 A dc from graph P.5.1, I, =53 600 A dc and, analytically, dc R d c 1.8670 J 3 t ) / A > Open circuit simulations: With the receiving end of the line open and i n i t i a l conditions equal to zero, a 60-Hz, l-p.u.-peak, sinusoidal voltage source was applied to the sending end of the line (peak voltage at t = 0). The resulting voltages at the receiving end are shown in the following graphs: i) Frequency dependence vs. 60 Hz •+ P.7.11 i i ) Frequency dependence vs. 750 Hz -»• P.7.12 For these simulations At = 0.01 msec. From the graphs, the maximum transient voltages are: i) Frequency-dependence model: 1.82 p.u. i i ) 60-Hz constant-parameter model: 1.95 p.u. (+7%) i i i ) 750-Hz constant-parameter model: 1.55 p.u. (-15%) The difference in the harmonic content is quite noticeable. The 60-Hz model exaggerates the harmonics above 60 Hz, and the 750-Hz model, while exaggerating the higher frequency harmonics, overdampens the mid-range ones. In a second open circuit simulation, the line was energized with a unit step voltage. The voltages at the receiving end are shown in the graphs: i) Frequency dependence vs. 60 Hz •+ P.7.13 i i ) Frequency dependence vs. 750 Hz P.7.14 (At = 0.01 msec). The peak voltages for this case were: 118 i ) Frequency dependence: 1.93 p.u. i i ) 60-Hz parameters: 1.96 p.u. (+1.6%) i i i ) 750-Hz pa r a m e t e r s : 1.60 p.u. (-17%) The 60-Hz model reproduced the maximum peak v o l t a g e q u i t e c l o s e l y , but t h e damping of the o s c i l l a t i o n s i s v e r y low. The 750-Hz model g i v e s an o v e r a l l e x c e s s i v e damping, even though i t e x a g g e r a t e s the v e r y sharp c o r n e r s ( h i g h e s t f r e q u e n c i e s ) . A d d i t i o n a l examples: As a d d i t i o n a l i l l u s t r a t i o n s o f the importance o f the c o r r e c t s i m u l a -t i o n o f the d i f f e r e n t f r e q u e n c y components of the s i g n a l i n t r a n s i e n t s i m u l a -t i o n s , the f o l l o w i n g examples a r e p r e s e n t e d . A) L i n e t e r m i n a t e d i n a c a p a c i t a n c e : In these s i m u l a t i o n s t h e l i n e was t e r m i n a t e d i n a c a p a c i t i v e l o a d e q u a l t o ^ the c a p a c i t a n c e of t h e l i n e , C l o a d = 7- C l i n e = (0.01244 — ) x 100 mi = 0.311 uF 4 4 mi 1 3 5 With At = 0.01 msec, to d i s c r i m i n a t e f r e q u e n c i e s up t o f = n n ^ x 10 = 1 0 Hz, the f o l l o w i n g s i m u l a t i o n s were r u n : a) For a u n i t s t e p - v o l t a g e i n p u t and z e r o i n i t i a l c o n d i t i o n s , t h e v o l t a g e a c r o s s t h e c a p a c i t a n c e i s shown i n the f o l l o w i n g g r a p h s : i ) Frequency dependence v s . 60 Hz -*• P.7.15 i i ) Frequency dependence v s . 750 Hz -»• P.7.16 b) For a s t e p v o l t a g e of 100 kV, t h e c u r r e n t i n t h e c a p a c i t a n c e i s shown i n t h e f o l l o w i n g g r a p h s : i ) Frequency dependence v s . 60 Hz ->- P.7.17 i i ) Frequency dependence v s . 750 Hz -*• P.7.18 119 B) Line terminated in an L-C load: C load = 0.311 uF L load = 13 mH The value of L load was selected to give a resonant frequency in the LC combination equal to 2500 Hz. For a step voltage of 100 kV applied at the beginning of the line, the following simulations were obtained (At = 0.01 msec): a) Load voltages (in p.u.): i) Frequency dependence vs. 60 Hz ->- P.7.19 i i ) Frequency dependence vs. 750 Hz P.7.20 b) Currents in the capacitance (in amps): i) Frequency dependence vs. 60 Hz ->• P.7.21 ii ) Frequency dependence vs. 750 Hz -> P.7.22 c) Currents in the inductance (in amps) : i) Frequency dependence vs. 60 Hz -* P.7.23 i i ) Frequency dependence vs. 750 Hz -* P.7.24 6.8 Field Test Simulation Oscillographs and data of a f i e l d test have been made available by the Bonneville Power Administration (BPA), Portland, Oregon, U.S.A. The test consisted of a single-line-to-ground fault on an open-ended, 222 km, 500 kV, 3-phase transmission lin e . The short circuit was applied to phase-c. The oscillographs of the voltages at the fault location are shown in graph P.8.1(a) . 120 The simulation of t h i s f i e l d test using BPA's weighting functions formulation i s discussed i n reference [13]. In t h i s simulation a s i m p l i f i e d representation was used for the generators and transformers feeding the l i n e This representation i s shown in f i g . 6.2. The tower footing resistance was assumed to be 2 fi and the t y p i c a l value of 100 fi-m was assumed for the earth r e s i s t i v i t y . 500 kV. 222 km F i g . 6.2: John Day-Lower Monumental BPA's f i e l d t e s t . Ref. [13] representation. In BPA's simulation the zero sequence mode of the l i n e was modelled for frequency dependence. The i n t e g r a t i o n step for the EMTP was At = 50 usee. To compare the r e s u l t s obtained with the l i n e model developed i n t h i s work with the r e s u l t s obtained in BPA's simulation, the same conditions as i n BPA's simulation were assumed. The phase voltages obtained with the new model are shown i n graph P.8.1(b). The voltages at phase-b corresponding to the f i e l d o s c i l l o g r a p h , BPA's simulation, and new model simulation are shown i n graphs P.8.2. As i t can be seen from graph P.8.2, there i s a good agreement between the simulations and the f i e l d o s c i l l o g r a p h . The peak voltages ( i n per u n i t of the p r e - f a u l t voltage at the s i t e of the f a u l t ) for the d i f f e r e n t cases (including a constant-60-Hz-parameter simulation) were: 121 i) Oscillograph: 1.60 p.u. i i ) BPA's simulation: 1.77 p.u. (+11%) i i i ) New model simulation: 1.71 p.u. (+7%) iv) 60-Hz constant parameters: 2.11 p.u. (+32%) The slightly higher value of the peak voltage given by the new model simulation as compared with the f i e l d test result may be due to a number of reasons, as for instance: a) A very simple source representation . b) Earth r e s i s t i v i t y different from assumed value. c) Frequency dependence of the positive sequence mode. These factors could also account for the slightly higher harmonic content of the frequency-dependence simulations as compared with the f i e l d test oscillograph. An interesting factor that can be observed in graphs P.8.2 is that in the case of the f i e l d test and new model simulation, the peak voltages tend to decrease as the steady state is approached. In BPA's simulation, however, these values seem to follow a slight increase. In BPA's frequency-dependence simulation the average time per step of the solution process was 3.13 times longer than with a constant-parameter representation. In the simulation with the new model this time was only 1.19 times longer. 122 (CHAPTER 6: NUMERICAL RESULTS) LIST OF GRAPHS LINE PARAMETERS P.2.1: Frequency dependence of the series resistance 127 P.2.2: Frequency dependence of the series inductance 127 CHARACTERISTIC IMPEDANCE Zero sequence: i) Form of Z c(OJ) P.3.1: Magnitude of Zc(w). Zero sequence 128 P.3.2: Phase angle of Z C ( O J ) . Zero sequence 128 P.3.3: Real part of Zc(co) . Zero sequence 12 9 P.3.4: Imaginary part of Z c ( u j ) . Zero sequence 129 i i ) Approximation of Z c ( c o ) P.3.5: Approximation of Z c(io), zero seq. Magnitude 130 P. 3.6: Approximation of Z C (OJ), zero seq. Phase angle 130 P.3.7: Approximation of Zc(w), zero seq. Real part 131 P.3.8: Approximation of Zc(co), zero zeq. Imaginary part 131 P.3.9: Approximation of Zc(w), zero seq. Magnitude error 132 P.3.10: Approximation of Zc(tu), zero seq. Phase error 132 Positive Sequence: i) Form of Zc(w) P.3.11: Magnitude of Zc(u>) . Positive sequence 133 P.3.12: Phase angle of Z c(u). Positive sequence 133 P.3.13: Real part of Zc(co) . Positive sequence 134 P.3.14: Imaginary part of Z c(OJ) . Positive sequence 134 i i ) Approximation of Z c(OJ) P.3.15: Approximation of Z c(OJ) , pos. seq. Magnitude 135 P.3.16: Approximation of Z c(OJ), pos. seq. Phase angle 135 P.3.17: Approximation of Z c(OJ), pos. seq. Real part 136 P.3.18: Approximation of Zc(w), pos. seq. Imaginary part 136 P.3.19: Approximation of Zc(w), pos. seq. Magnitude error 137 P.3.20: Approximation of Zc(w), pos. seq. Phase error 137 WEIGHTING FUNCTION Simulation of A^(w) for Different Line Lengths P.4.1: Simulation of A-j_(co) . Zero sequence, 5-mi 138 P.4.2: Simulation of A^ (oo) . Positive sequence, 5-mi 138 P.4.3: Simulation of A^(OJ) . Zero sequence, 30-mi 139 123 P.4.4: Simulation of A-j_(io). Positive sequence, 30-mi 139 P.4.5: Simulation of A-i (to) . Zero sequence, 100-mi 140 P.4.6: Simulation of A-L(CO) . Positive sequence, 100-mi 140 P.4.7: Simulation of A-j_(to) . Zero sequence, 500-mi 141 P.4.8: Simulation of A^(co) . Positive sequence, 500-mi 141 P.4.9: Simulation of A-j_(ti)) . Zero sequence, 1000-mi 142 P.4.10: Simulation of A-^Cto) . Positive sequence, 1000-mi 142 Zero Sequence, 100-mi: i) Form of A-i_(co) P.4.11: Magnitude of A^(co) . Zero sequence, 100-mi 143 P.4.12: Phase angle of A]_(a)) . Zero sequence, 100-mi 143 P.4.13: Real part of A-J_(LO) . Zero sequence, 100-mi 144 P.4.14: Imaginary part of A^Cto) . Zero sequence, 100-mi 144 P.4.15: Phase angle of P(co) . Zero sequence, 100-mi 145 i i ) Simulation of A-L(IJJ) P.4.16: Simulation of A-j_(cu), zero seq., 100-mi. Magnitude 146 P.4.17: Simulation of A^(w) , zero seq., 100-mi. Phase angle 14 6 P.4.18: Simulation of A;j_(tu), zero seq., 100-mi. Real part 147 P.4.19: Simulation of A-j^ (oa) , zero seq., 100-mi. Imaginary part 147 P.4.20: Simulation of A^Cto) , zero seq., 100-mi. Magnitude error 148 P.4.21: Simulation of A-]_(a)) , zero seq., 100-mi. Phase error 148 Positive Sequence, 100-mi: i) Form of A^(cu) P.4.22: Magnitude of A^Cto) . Positive sequence, 100-mi 149 P.4.23: Phase angle of A]_(CJ) . Positive sequence, 100-mi 149 P.4.24: Real part of A^(io) . Positive sequence, 100-mi 150 P.4.25: Imaginary part of A^Cco) . Positive sequence, 100-mi 150 P.4.26: Phase angle of P(co) . Positive sequence, 100-mi 151 i i ) Simulation of A^(w) P. 4.27: Simulation of A^(tu), pos. seq., 100-mi. Magnitude 152 P.4.28: Simulation of A-j_(co) , pos. seq., 100-mi. Phase angle 152 P.4.29: Simulation of A-]_((JJ) , pos. seq., 100-mi. Real part 153 P.4.30: Simulation of A-j_(co), pos. seq., 100-mi. Imaginary part 153 P.4.31: Simulation of A^(u}) , pos. seq., 100-mi. Magnitude error 154 P.4.32: Simulation of A^(a)) , pos. seq., 100-mi. Phase error 154 Time Domain Form, a ^ C t ) : P.4.33: Weighting function a^(t). Zero seq., 100-mi 155 P.4.34(a):Weighting function a-]_(t). Pos. seq., 100-mi. 155 Same time scale as in graph P.4.33. P. 4. 34(b): Weight ing function a-j_(t). Pos. seq., 100-mi. 155 Expanded time scale. 124 FREQUENCY DOMAIN TESTS Zero Sequence: i ) S h o r t - c i r c u i t r e s p o n s e s p . 5 . 1 ( a ) : S/C Response, z e r o seq. h i g h f r e q u e n c i e s . E x a c t v s . approx. Mid to 156 p .5.1(b) : S/C Response, z e r o seq. f r e q u e n c i e s . E x a c t v s . approx. Low 156 p . 5 . 2 ( a ) : S/C Response, z e r o seq. to h i g h f r e q u e n c i e s . E x a c t v s . 60-Hz. Mid 157 p .5.2(b) : S/C Response, z e r o seq. f r e q u e n c i e s . E x a c t v s . 60-Hz. Low 157 p .5.3(a) : S/C Response, z e r o seq. h i g h f r e q u e n c i e s . E x a c t v s . 750-Hz. Mid t o 158 p .5.3(b) : S/C Response, z e r o seq. E x a c t v s . 750-Hz. Low 158 .5.4(a) : f r e q u e n c i e s . p S/C Response, z e r o seq. Magnitude e r r o r . A p p r o x i m a t i o n . 159 Mid to h i g h f r e q u e n c i e s . p .5.5(a) : S/C Response, z e r o seq. Magnitude e r r o r . 60-Hz. 159 .5.6(a) : Mid to h i g h f r e q u e n c i e s . p S/C Response, z e r o seq. Mid to h i g h f r e q u e n c i e s . Magnitude e r r o r . 750-Hz. 159 p .5.4(b) : S/C Response, z e r o seq. Low f r e q u e n c i e s . Magnitude e r r o r . Approx imat i o n . 160 p .5.5(b) : S/C Response, z e r o seq. Low f r e q u e n c i e s . Magnitude e r r o r . 60 Hz. 160 p .5.6(b) : S/C Response, zero seq. Low f r e q u e n c i e s . Magnitude e r r o r . 750 Hz. 160 i i ) O p e n - c i r c u i t r e s p o n s e s P.5. 7 : 0/C Response, zero seq. Exact vs. approx. 161 P.5. 8: 0/C Response, zero seq. Exact vs. 60-Hz. 161 P.5. 9: 0/C Response, zero seq. Exact vs. 750-Hz. 162 P.5. 10: 0/C Response, zero seq. Magnitude error. Approximat ion. 163 P.5. 11: 0/C Response, zero seq. Magnitude error. 60-Hz. 163 P.5. 12: 0/C Response, zero seq. Magnitude error. 750-Hz. 163 P o s i t i v e Sequence: i ) S h o r t - c i r c u i t r e s p o n s e s P.5.13(a):S/C Response, pos. seq. E x a c t v s . approx. Mid to 164 h i g h f r e q u e n c i e s . P.5.13(b) :S/C Response, pos. f r e q u e n c i e s . seq. E x a c t v s . approx. Low 164 P.5.14(a) :S/C Response, pos. h i g h f r e q u e n c i e s . seq. E x a c t v s . approx. Mid to 165 P.5.14(b) :S/C Response, pos. f r e q u e n c i e s . seq. E x a c t v s . 60-Hz. Low 165 P.5.15(a) :S/C Response, pos. h i g h f r e q u e n c i e s . seq. E x a c t v s . 750-Hz. Mid to 166 P.5.15(b) :S/C Response, pos. f r e q u e n c i e s . seq. E x a c t v s . 750-Hz. Low 166 P.5.16(a) :S/C Response, pos. seq. Magnitude e r r o r . A p p r o x i m a t i o n . 167 Mid t o h i g h f r e q u e n c i e s . 125 p •5.17(a) :S/C Response, pos, seq. Magnitude error. 60-Hz. 167 Mid to high frequencies. p .5.18(a) :S/C Response, pos. seq. Magnitude error. 750-Hz. 167 Mid to high frequencies. p .5.16(b) :S/C Response, pos. seq. Magnitude error. Approximation. 168 Low frequencies. p .5.17(b) :S/C Response, pos. seq. Magnitude error. 60-Hz. 168 Low frequencies. p .5.18(b) :S/C Response, pos. seq. Magnitude error. 750-Hz. 168 Low frequencies. Open-circuit responses P .5.19 0/C Response, pos. seq. Exact vs. approx. 169 P .5.20 0/C Response, pos. seq. Exact vs. 60-Hz. 169 P .5.21 O/C Response, pos. seq. Exact vs. 750-Hz. 170 P .5.22 0/C Response, pos. seq. Magnitude error. Approximat ion. 171 P .5.23 0/C Response, pos. seq. Magnitude error. 60-Hz. 171 P .5.24 0/C Response, pos. seq. Magnitude error. 750-Hz. 171 TRANSIENT SIMULATIONS Short Circuit: P. 7 .1: Short circuit simulation. Sinusoidal source , zero at 172 t=0. New model. p. 7 .2: Short circuit simulation. Sinusoidal source , zero at 172 t=0. 60-Hz parameters. p. 7 .3 : Short circuit simulation. Sinusoidal source , zero at 173 t=0. 750-Hz parameters. p. 7 .4: Short circuit simulation. Sinusoidal source , zero at 174 t=0. New model vs. 60-Hz. First cycles. p. 7 .5: Short circuit simulation. Sinusoidal source , zero at 174 t=0. New model vs. 750-Hz First cycles. p. 7 .6: Short circuit simulation. Sinusoidal source , peak at 175 t=0. New model vs. 60-Hz. p. 7 .7: Short circuit simulation. Sinusoidal source , peak at 175 t=0. New model vs. 750-Hz p. 7 .8: Short circuit simulation. Step excitation. New model 176 vs. 60-Hz. p. 7 .9: Short circuit simulation. Step excitation. New model 17 6 vs. 750-Hz. p. 7 .10: Short circuit simulation. Step excitation. New model. 177 Extended time response. Open Circuit: P.7.11: Open cir c u i t simulation. Sinusoidal source, peak at 178 t=0. New model vs. 60-Hz. P.7.12: Open circuit simulation. Sinusoidal source, peak at 178 t=0. New model vs. 750 Hz. P.7.13: Open cir c u i t simulation. Step excitation. New model 179 vs. 60-Hz. 12 6 P.7.14: Open circuit simulation. Step excitation. New model 179 vs. 750-Hz. Capacitive Loads: P. 7 .15: Capacitive load. Step New model vs. 60-Hz. excitation. Load voltage. 180 p. 7 .16: Capacitive load. Step New model vs. 750-Hz. excitation. Load voltage. 180 p. 7 .17: Capacitive load. Step New model vs. 60-Hz. excitation. Load current. 181 p. 7 .18: Capacitive load. Step New model vs. 750-Hz. excitation. Load current. 181 L-C Load: P. 7 .19: L-C load. Step excitation. Load voltage. New 182 model vs. 60-Hz. p. 7 .20: L-C load. Step excitat ion. Load voltage. New 182 model vs. 750-Hz p. 7 .21: L-C load. Step excitation. Current in C. New 183 model vs. 60-Hz. p. 7 .22: L-C load. Step excitation. Current in C. New 183 model vs. 750-Hz p. 7 .23 : L-C load. Step excitation. Current in L. New 184 model vs. 60-Hz. p. 7 .24: L-C load. Step excitation. Current in L. New 184 model vs. 750-Hz m FIELD TEST SIMULATIONS P.8.1(a): Field test oscillograph (BPA). Single-line-to-ground 185 fault on phase C. Receiving end voltages. P.8.1(b): Simulation of f i e l d test in P.8.1(a) with the new line 185 model. (Only the zero sequence is modelled as frequency dependent.) P.8.2: Comparison between voltages at phase b for: a) Field 186 test oscillograph; b) BPA's frequency-dependence simulation; c) New model simulation. 127 1(T q IO 3 A g]0 2 • \ to T. 210 UJ <_> 1 d CC I — t o to LU CE]0-» 1 0 2 4 . i i in II—i i in II—i i in II—i n u n — i i in II—i i in II—i i in II—i n u n — i i in II—i i in II—i i in II IO"' 10"2 10"' ] 10 10* 10' 10' 10s 106 JO7 10" FREOUENCT (HZ) P.2.1: Frequency dependence of the series resistance. 8 - j 7 -= 6 H io <• u -z. CC o 3 H ZD a -z. ~ 2 1 Pos. S e q . i i III i i — i i III II—i i III i i — i mm—i i in II—i i in II—i i in II—i mm—i i in n — i i IIMI—rrrrrn 10"3 10"? 10"' 1 10 IO1 103 1 0' 10s 106 1 07 !0 8 FREQUENCY (HZ) P.2.2: Frequency dependence of the series inductance. 128 P.3.2: Phase angle of Z (to). Zero sequence. -80 H •100 I i i i i in IIIIIII i imn i 11 mi l 111111 i min i nun i min i nun i nun i nun _ 10"* 10"* 10"' 1 10 10* 1 Cf 10' 10s 106 10 10e FREQUENCr (HZ) P.3.4: Imaginary part of Z (to). Zero sequence. 130 P. 3.6: Approximation of Z (to), zero seq. Phase angle. 131 800 -| P.3.7: A p p r o x i m a t i o n of Z (to), zero seq. R e a l p a r t . P.3.8: A p p r o x i m a t i o n o f Z (to), z e r o s e q . Imaginary p a r t . 0 I I I I N I I i i i n i i 11 I I I ii i i iM u i T 11 T 11 i m i n i nun i nun i m i i i— i uu u i i in 11 10 i c f 10"' i io io ! icf io" io 5 io' io 7 icf FREQUENCY (HZ) 9: Approximation of Zc(co), zero seq. Magnitude error. 0 0 (H 0 0 0 0 1 I I I I I I I I T i II i < i nun i i in II i lui I I i n u n i i i n i i i H U M I I I I I M I I I I I I I I I I I I I I 10 10"' 10"' 1 10 10' l l f 10" 10s 10" io 1 io" FREOUENCT (HZ) 3.10: Approximation of Z^(ai), zero seq. Phase error. 133 P.3.12: Phase angle of Z (to) . Positive sequence. 200 I IIIIIII l imn l imn i mm M i n n i i in n l imn i nun l imn l imn M i n n 10"3 10 10"' 1 10 102 IO3 10' 10 10 10 10 FREQUENCY (HZ) P.3.14: Imaginary part of Z ( t o ) . Positive sequence. 135 700 to z: x o o LU CO to O Q_ O U J to M 300 A 250 -1 10' IIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII i i in II i i in n 10" 10"' 1 10 i o 2 l ( f 10' 10 10 10 10 FREOUENCT (HZ) P.3.15: Approximation of Z^Cco), pos. seq. Magnitude. LO o o LU t o to a o UJ CL O <_) X a. 30 25 -20 -15 -j 10 -j 5 - i 0 -5 --]0 ~ -15 -20 --25 --30 10' IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII I IIIIII I IIIIII I I III II I I III II 10' 10 10 10 10 10 FREOUENCT (HZ) 10T 10 10 10 P.3.16: Approximation of Z^(ud) , pos. seq. Phase angle. 11 III II II III II i i III II i min i i in II—i i in II—i i in II II in II—i i in II—i i in II i i in n 10" 10 10 1 10 10 i o 3 10' i o 5 10* i o 7 i o 8 FREQUENCY (HZ) P.3.17: Approximation of Z c(OJ), pos. seq. Real part. -200 i i III II 11 in II j i in II—i min—i i in II—i i in II—i i in II i i nm i i m n—i mm—11 in ri 10 10 10 1 10 10 103 10' 10s 106 101 10e FREQUENCY (HZ) P.3.18: Approximation of Z (co), pos. seq. Imaginary part. 137 5.0 i 4.0 3.0 2.0 1.0-J 0 -1 .0 -2.0 -3.0 : -4.0 --5.0 IIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII IIIIIII I I III II I I III II 10"3 10 10 1 10 10* 103 10" 10 10 10 1 or FREOUENCT (HZ) P.3.19: Approximation of Z c(to), pos. seq. Magnitude error, 5 0 o o o c 4 0 -o or : CC U J 3 0 -o L U 2 0 -t o U-> 1 0 -O CL ,—, 0 a ru -1 0 -X CL -2 0 -CO > , -3 0 -o X -4 0 -a_ -5 0 I I III II II III II I I III II I I Mill I I III II—I I III 11—I IIIIII I IIIIII I IIIIII I IIIIII I IIIIII 10 10 10"' 1 10 10! 103 10" 10s 10* io 1 icf FREOUENCT (HZ) P.3.20: Approximation of Z (to), pos. seq. Phase error. 138 1 .0 0 • 9 : o LJ 0 • 8 : to • o oc 0 7 : LU rsl 0 6 : Ti 1 1 in 0 5 : 1 Pll 0 4 \ VS. 0 3 : lib 1 0 2 : 0 1 : 0.0 11 III II II IM II i i III II—i n u n — i i in II—i i in II—i i in II—i HUM—i i MI n r i III II—i i m n 10 10"2 10"' 1 10 102 l l f 10' 1 05 1 06 1 07 1 0e FREQUENCY (HZ) P.4.1: S i m u l a t i o n o f A^w) . Zero sequence, 5-mi. A1 O UJ t o to o a. IT) to > c r 1 .0 0.9 0.8 0.7H! 0.6 0.5 -| 0.4 0.3 -0.2 : 0.1 : 0 • 0 I II III II II III II I i III II i mm—i i in n — i i in n — i i in II—i mm—i i in n — i MINI—i i in n 10 10 10 1 10 10 103 10' 1 05 1 06 1 07 1 0" FREQUENCY (HZ) P.4.2: S i m u l a t i o n o f A ^ t o ) . P o s i t i v e sequence, 5-mi. 139 FREOUENCT (HZ) P.4.3: Simulation of A., (to) . Zero sequence, 30-mi. FREOUENCT (HZ) P.4.4: Sim u l a t i o n of A., (to) . P o s i t i v e sequence, 30-mi. 140 P.4.6: S i m u l a t i o n of A., (co) . P o s i t i v e sequence, 100-mi. FREQUENCY (HZ) P.4.7: S i m u l a t i o n o f A., (to) . Zero sequence, 500-mi. A l 1 . 0 ^ FREQUENCY (HZ) P.4.8: S i m u l a t i o n of A (to) . P o s i t i v e sequence, 500-mi. 142 FREQUENCY (HZ) P.4.9: S i m u l a t i o n of A., (to) . Zero sequence, 1000-mi. FREQUENCY (HZ) P.4.10: S i m u l a t i o n of A., (to) . P o s i t i v e sequence, 1000-mi. 143 FREQUENCY (HZ) P.4.11: Mangitude of A., (OJ) . Zero sequence, 100-mi. FREQUENCY (HZ) P.4.12: Phase angle of A (to) . Zero sequence, 100-mi. 144 i I I I I I I — i — i i n i i 1—i III II—i—i M I II 1—i i n II 1—IIIIII 1 0 2 4 1 2 4 1 0 2 4 1 0 2 4 1 o' 2 4 10" 2 4 1 0 S FREOUENCT (HZ) P.4.13: Real part of A^to) . Zero sequence, 100-mi. ~i—i III II 1—i III l i 1—i i n l l 1—IIIIII 1—I i n l i 1—i i n II 10" 1 2 4 1 2 4 1 0 2 4 1 0 2 2 4 1 0 S 2 4 10* 2 4 1 0 S FREOUENCT (HZ) P.4.14: Imaginary part of A^(co) . Zero sequence, 100-mi. 145 P.4.15: Phase angle of P(co) . Zero sequence, 100-mi. .146 o UJ CO o m LU M I o o CO > cr 1 .0 0.9 0.8 0.7 0.6 : 0.5 : 0.4 : 0.3 : 0.2 : 0 . 1 -"0 .0 I i i III II i i III II i i III II IIIIIII I I I ,-3 , ^ -2 10 10' 10 i i IM II—i i ill n r i in I I—i i HI I I—i i in I I—i i in n 10 10 10 10" 10s 106 101 108 FREQUENCY (HZ) P.4.16: S i m u l a t i o n of A ^ t o ) , z e r o seq., 100-mi. Magnitude. o CC UJ rsl X 0_ CO cr x Q . T 1 I I I I I I I 1 1 I I I I III 1 1 IIIIIII 2 3 4 6 10 2 3 4 6 10* 2 3 4 6 103 2 3 4 6 10" 2 3 4 6 10s FREQUENCY (HZ) P.4.17: S i m u l a t i o n o f A^ Cw), z e r o s e q . , 100-mi. Phase a n g l e . 147 1 .0 7 0.8 : 0.6 : 0.4 : 0.2 : •0.0 : •0.2 -0.4 -•0.6 --0.8 : 1 .0 H 1—I I I I II 1—I M l II 1—I III II 1—I III II 1—I 11 I 11 1—I III II 10"'2 4 1 2 4 10 2 4 ]02 2 4 10' 2 4 10* 2 4 10s FREQUENCY (HZ) P.4.18: Simulation of A 1 ( u ) , zero seq., 100-mi. Real part. 1 .0 -; 0.8 \ 0.6 \ 0.4^ 0.2 : 0 .0 -•0.2 -•0.4 -j •0.6 \ •0.8 \ •1.0 -1 1 I I I I I I 1 I I I I I I 1 I I I I I I i I I I I I I , i I I I I I I i I I I I I I 10"'2 4 1 2 4 10 2 4 10 2 4 10 2 4 10 2 4 10 FREQUENCY (HZ) P.4.19: Simulation of A (to), zero seq., 100-mi. Imaginary part. 148 P.4.20: Simulation of A (ID) , zero seq., 100-mi. Magnitude error. ID Q CC o ce cc UJ o UJ • 0 i .0 .0 .0 • CH 0 ± -1-X Q. • -2 i n > - -3 cr i -A CH 0 0 : 0 -0 I I I I I I — I I I I I I I — I I I I I I I — I I I I I I I — I I I I I I I — I I I I I I I — I I I I I I I — I I I I I I I 1(T 10" 10"' I 10 icf icf io" icf FREQUENCY (HZ) P.4.21: Simulation of A (ID) , zero seq., 100-mi. Phase error. 149 FREQUENCY (HZ) P.4.23: Phase angle of AALO) . Positive sequence, 100-mi. 150 1 .0 10"' 1 10 10* 10' 10' 10s 10* 1 07 3 0° FREQUENCY (HZ) P.4.24: R e a l p a r t of A 1 (to) . P o s i t i v e sequence, 100-mi. FREQUENCY (HZ) P.4.25: Imaginary p a r t o f A (to) . P o s i t i v e sequence, 100-mi. 151 P.A.26: Phase angle of P(OJ) . Positive sequence, 100-mi. A1 P.4.28: Simulation of A (to) , pos. seq., 100-mi. Phase angle. 153 FREQUENCY (HZ) P.4.29: Simulation of A (to) , pos. seq., 100-mi. Real part. FREQUENCY (HZ) P.4 .30: Simulation of A., (to) , pos. seq., 100-mi. Imaginary part. 154 oc o o o Q. 5.0 4.0 3.0 2.0 1.0 d 0 -1.0 -J -2.0 -3.0 -4.0 -1 -5.0 I Ml II J I III II—I I III II—I l l l l II—I I 11I 11—I I III II—I I III II—I 10"J 10 10" 1 10 101 ltf 10' FREQUENCY (HZ) 10' r r m — i I m n 10' 101 P.4.31: Simulation of A^ (oj) , pos. seq., 100-mi. Magnitude error. o o a: a cn CK U J x Q_ 5.0 : 4.0 -3.0 : 2.0 -1.0 -0 -1.0 -2.0 -3.0 -4.0 -I -5.0 10' Mini—i i M I I I—i i in I I 10"' 10"' 1 TTT71 1 M U M 1 I III II 1 I III II 1 M i l II 1 I 111 TI 1 I III II 10 10! 10s 10' 10s io" io 7 FREQUENCY (HZ) P.4.32: Simulation of A (LO) , pos. seq., 100-mi. Phase error. 6000 o 5000 UJ in co ce & 4000 5 3000 13 2 2000 x L3 1000 0.5 1.0 1.5 2.0 TIME (MSEC) 2.5 3.0 F.4.33: Weighting function a ^ t ) . Zero seq., 100-mi. o u i V) a a. i o U J 400000 350000 300000 : 250000 -200000 150000 100000 50000 H" 0.5 1 .0 1.5 2.0 TIME (MSEC) 2.5 3.0 P.4.34(a): Weighting function a j ^ t ) . Pos. seq., 100-mi. Same time scale as i n graph P.4.33. o U J co o o_ o I— CJ z U J 3 400000 350000 300000 250000 200000 150000 100000 50000 • • • i • ' 1 • i • • ' • i • • • • i • i • i | i i i • | • • i • | i • . . | 50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.60 TIME (MSEC) P.4.34(b): Weighting function a^(t) . Pos. seq., 100-mi. Expanded time scale. 156 800 1 700 A 5 600 A co > CJ o UJ CO 500 -] 400 A § 300 UJ 200 100 A 10 i i i i 111 1—i—i i i i 111 1—i—i II i i I I — 2 3 4 6 103 2 3 4 6 10* FREQUENCY (HZ) 2 3 4 6 10 -i 1—i i i i II i 2 3 4 6 10s P.5.1(a): S/C Response, zero seq. Exact vs. approx. Mid to high frequencies. 60000 -| FREQUENCY (HZ) P.5.1(b): S/C Response, zero seq. Exact vs. approx. Low frequencies. 157 5000 n 5 4000 H o t£) 0. 2 3000 -o FREQUENCT (HZ) P.5.2(a): S/C Response, zero seq. Exact vs. 60-Hz. Mid to high frequencies. P.5.2(b): S/C Response, zero seq. Exact vs. 60-1 Low frequencies. 158 800 -) 0 H. 1 1 1 I I I I I I 1 1 1 I I I I I I 1 1 1 I M i l l 1 1 1 I I I I I I 10 2 3 4 6 10J 2 3 4 6 103 2 3 4 6 10* 2 3 4 6 10 FREOUENCT (HZ) " P . 5 . 3 ( a ) : S/C Response, z e r o seq. E x a c t v s . 750-Hz. Mid to h i g h f r e q u e n c i e s . 60000 n FREQUENCY IHZI P.5.3(b): S/C Response, z e r o seq. E x a c t v s . 750-Hz. Low f r e q u e n c i e s . at o oc 5 in > o CE 10 1 i I I I I I I I 1 i I I I I I I I 1 i I I I I I I I 1 — i — i i i 1111 2 3 4 6 itf 2 3 4 6 itf 2 3 4 6 10' 2 3 4 6 JO* FREOUENCT (HZ) P.5.4(a): S/C Response, zero seq. Magnitude error. Approximation. Mid to high frequencies. 500 -400 -•* OC 300 : o 1 200 -o IO 100 -a. in > o -u -100 -a UJ in -200 -o oc UJ M -300 : 3C -400 --500 : -I 1 I I I I I I I . 1 1 I I I Mil 1 1 I I I I I I I 1 1 I I I I I I I 10 2 3 4 6 ltf 2 3 4 6 i t f - 2 3 4 6 10 2 3 4 6 10* FREOUENCT (HZ) P.5.5(a) : S/C Response, zero seq. Magnitude error. 60-Hz. Mid to high frequencies. oc D a 5? r-a. in o UJ in o oc U J M 10 I 1 I I I I I I I . 2 3 4 6 10* T 1 I II I I II 1 1 I I I I I I I 2 3 4 6 l(f 2 3 4 6 10* FREOUENCT (HZ) I I I i m i 2 3 4 6 10* .5.6(a): S/C Response, zero seq. Magnitude erro r . Mid to high frequencies. 750-Hz. 5 0 o ce 4 0 3 0 : 2 0 : 10 0 - 1 0 - 2 0 - 3 0 - 4 0 -50 —i— i I I I I I I I 1—i I I I I I I I 1—i I I I I I I I 1—i I I I I I I I 1—i i I I i in 10"'' 2 3 4 6 10"' 2 3 4 6 10"' 2 3 4 6 1 2 3 4 6 10 2 3 4 6 10* FREOUENCT (HZ) P.5.4(b): S/C Response, zero seq. Magnitude erro r . Approximation. Low frequencies. g 5 2 0. in >• 100 eo 60 40 20 -20 -40 -j -60 -80 -J -100 P60 10" 8 o. a UJ m o tr. 100 80 60 -40 20 -0 -20 -40 -60 -] -80 -100 . 1 4 Vc rft X 4 ! • " " " -i ' " ' i ' m m 1 i i m m 2 34 6 10 2 3 4 6 10' 2 34 6 1 2 3 4 6 10 2 34 6 10* FREOUENCT (HZ) J P.5.5(b): S/C Response, zero seq. Magnitude error. 60-Hz. Low frequencies. 10" I I II I mi , I I I I 11 III . , I I I m i l 1 l I I I m i . 2 34 6 10' 2 3 4 6 10' 2 34 6 1 2 3 4 6 10 2 3 4 6 10* FREOUENCT (HZ) .5.6(b): S/C Response, zero seq. Magnitude error. 750-Hz. Low frequencies. 161 10 1 FREOUENCT (HZ) P.5.7: 0/C Response, zero seq. Exact vs. approx. 40 n 35 : ZD 0 H 1 1 1 I I I I I I 1 1 1 I I I I I I 1 1 I I I M I T — I I I I I I I I I 10 2 3 4 6 102 2 3 4 6 103 2 3 4 6 10* 2 3 4 6 10 FREOUENCT (HZ) P.5.8: 0/C Response, zero seq. Exact vs. 60-Hz. 162 10 4 6 10 3 4 6 10 FREQUENCY 2 3 (HZ) P.5.9: 0/C Response, zero seq. Exact vs. 750-Hz. 10 P60 "10 1 i i — i i I I i II — i — i — i i i 1111 1 — i — i i i 1111 1 — i — i i i 1111 10 2 3 4 6 10 2 3 4 6 itf 2 3 4 6 10* 2 3 4 6 10s FREQUENCY !HZ) P.5.10: 0/C Response, zero seq. Magnitude error. Approximation. 500 400 300 200 100 0 -100 -200 -300 -•400 -500 - 1 — i — i i i 1111 1 — i — i i i 1111 1 — i — i i i 1111 1 — i — i i i 1111 2 3 4 6 10! 2 3 4 6 ID3 2 3 4 6 10* 2 3 4 6 10* FREQUENCY (HZ) 10 P.5.11: O/C Response, zero seq. Magnitude error. 60-Hz. o UJ in o 500 400 300 200 100 0 -100 -200 •] -300 -400 -I -500 10 - i — i I I I I I I I , , , , , i -2 3 4 6 10! 2 3 4 6 Id3 2 3 ^ 6 JO" FREQUENCY (HZ) 2 3 4 6 10 P.5.12: 0/C Response, zero seq. Magnitude error. 750-Hz. 164 10 2 3 4 6 10 2 3 46 10 2 3 4 6 3 0 2 3 4 6 10 2 3 4 6 10 FREQUENCY (HZ) P.5.13(a): S/C Response, pos. seq. Exact vs. approx. Mid to high frequencies. 70000 i P.5.13(b): S/C Response, pos. seq. Exact vs. approx. Low frequencies. 165 FREQUENCY (HZ) P.5.14(a): S/C Response, pos. seq. E x a c t v s . 60-Hz. Mid to h i g h f r e q u e n c i e s . 70000 -i 60000 50000 40000 30000 A 20000 A I O O O O A —I—I I I I I I I I 1—I I I I I I I I 1—I I I I I I I I 1—I I I I I I I I 1—I I IIIIII 10'' 2 3 4 6 10"z 2 3 4 6 10"' 2 3 4 6 1 2 3 4 6 10 2 3 4 6 102 FREQUENCY (HZ) P.5.14(b): S/C Response, pos. seq. Ex a c t v s . 60-Hz. Low f r e q u e n c i e s . 166 10000 i i 1—i IIIIIII 1—i IIIIIII I i i i i l m 1 — I I l i n n 1—I i M I n i 10 2 3 4 6 llf 2 3 4 6 103 2 3 4 6 10* 2 3 4 6 ) 05 2 3 4 6 10* FREQUENCY (HZ) P.5.15(a): S/C Response, pos. seq. Exact vs. 750-Hz. Mid to high frequencies. 70000 -i 60000 A 50000 A 40000 A 30000 A 20000 A 10000 — I — I IIIIIII 10"3 2 3 4 6 10' -]—i IIIIIII 1—i IIIIIII 1—i IIIIIII 1—i IIIIIII 2 3 4 6 )0"' 2 3 4 6 1 2 3 4 6 10 2 3 4 6 1 02 FREQUENCY (HZ) P.5.15(b): S/C Response, pos. seq. Exact vs. 750-Hz. Low frequencies. 50 167 40 30 20 10 0 -10 -20 --30 --40 -50 10 P - i — i I I I I I I I - i — i I I I I I I I i— I , i i 11 , • •>' > i. i i i ini - i i i I m i l 1—i i i 11 i n 2 3 4 6 10 2 3 4 6 Jjf 2 3 4 6 JO 2 34 6 10 2 3 4 6 10s FREOUENCT (HZ) 5.16(a): S/C Response, pos. seq. Magnitude error. Approximation. Mid to high frequencies. CO ~i—i I I I I I I I 1—i I I I I I I I 1—i I I I I I I I 1—i I I I I I I I — — i — i I I I I I I I 10 2 3 4 6 I if 2 3 4 6 itf 2 3 4 6 IO* 2 3 4 6 10s 2 3 4 6 10* FREOUENCT (HZ) P.5.17(a): S/C Response, pos. seq. Magnitude error. 60-Hz. Mid to high frequencies. o U J to " i — I I I 11 in — I — I I I I I I I I — I — I i I I i m 1—i I I I I I I I 1 — i i i 11 in 10 2 34 6 1(T 2 3 4 6 10 2 3 4 6 10" 2 3 4 6 10s 2 34 6 io" FREOUENCT (HZ) P.5.18(a): S/C Response, pos. seq. Magnitude error. 750-Hz. Mid to high frequencies. o ce ce 100 80 60 40 20 0 -20-| -40 -60 -80 -!00 - 1 — i I I I I I I I — 10" 2 3 4 6 10"! P.5.16(b) ce ce UJ o UJ t o 100 60 60 40 20 0 -20 H -40 -60 -80 -j -100 ~ i — i i M m i 1—i I I I I I I I 1—i i i i n i l 1 — i I I I I I I I 2 3 4 6 10"' 2 3 4 6 1 2 3 4 6 10 2 3 4 6 10J FREOUENCT (HZ) S/C Response, pos. seq. Magnitude error. Approximation. Low frequencies. P60 — i — i I I I I I I I — !0" 2 3 4 6 10"! P.5.17(b): !00 80 60 40 20 0 -20 -40 -60 -60 -| 100 1 i i 111in 1—i I I I I I I I 1—i I I I I I I I 1—i I I I I I I I 2 3 4 6 10"' 2 3 4 6 1 2 3 4 6 !0 2 3 4 6 IO 2 FREOUENCT (HZ) S/C Response, pos. seq. Magnitude error. 60-Hz. Low frequencies. P750 - 1 — i I I I I I I I — 10"' 2 34 6 10"! P.5.18(b); i I I I I I I I 1 — r T l I I I I I I I I I 1—I I I I I I I I 2 3 4 6 10" 2 3 4 6 1 2 3 4 6 10 2 3 4 6 IO2 FREOUENCT (HZ! S/C Response, pos. seq. Magnitude error. 750-Hz. Low frequencies. 20 18 16 14 12 10 8 6 4 2 H — I — I IIIIIII 1—I IIIIIII 1—I IIIIIII 1—I IIIIIII 1—I IIIIIII 10 2 3 4 6 102 2 3 4 6 103 2 3 4 6 10* 2 3 4 6 10s 2 3 4 6 106 FREQUENCY (HZ) P.5.19: 0/C Response, pos. seq. Exact vs. approx. P.5.20: 0/C Response, pos. seq. Exact vs. 60-Hz. 170 P.5.21: 0/C Response, pos. seq. Exact vs. 750-Hz. 50 ce ce tu 40 30 -20 -10 : 0 -30 -20 -30 -40 -50 —I i I I I I I I I I I IIIIIII 1 1 IIIIIII 1 1 I I I MM 1 1 IIIIIII 10 2 3 4 6 IO2 2 3 4 6 10 2 3 4 6 IO" 2 3 4 6 io' 2 3 * 6 IO" FREQUENCY (HZ) P.5.22: 0/C Response, pos. seq. Magnitude erro r . Approximation. > . I IIIIIII 1 1 IIIIIII 1 1 IIIIIII 1 1 I I I MM 1 1 IIIIIII 10 2 3 4 6 10* 2 3 4 6 IO3 2 3 4 6 IO" 2 3 4 6 10s 2 3 4 6 IO* FREQUENCY (HZ) P.5.23: 0/C Response, pos. seq. Magnitude erro r . 60-Hz. ce o in o UJ to T 1 1 I I I I III 1 1 IIIIIII 1 1 IIIIIII 1 1 IIIIIII 1 1 I I I Mil 10 2 3 4 6 IO2 2 3 4 6 Jtf 2 3 4 6 10" 2 3 4 6 10s 2 3 4 6 IO* FRE0UENCY (HZ) P.5.24: O/C Response, pos. seq. Magnitude erro r . 750-Hz. 172 1000 n 600 O \ CO 200 -200 -600 H -1000 I • i 0 0.20 0.40 0.60 0.80 1.00 TIME (SEC) P.7.1: Short c i r c u i t simulation. Sinusoidal source, zero at t=0. New model. P.7.2: Short c i r c u i t simulation. Sinusoidal source, zero at t=0. 60-Hz parameters. 173 P.7.3: Short circuit simulation. Sinusoidal source zero at t=0. 750-Hz parameters. 174 1000 -] -600 i -800 --1000 1 i i -0.00 0.02 0.04 0.06 0.08 0.10 TIME (SEC) P.7.4: Short circuit simulation. Sinusoidal source, zero at t=0. New model vs. 60-Hz. First cycles. 1000 "] ^ -400 : -600 : -800 A -1000 I 1 1 ' i i i i -0 00 0.02 0.04 0.06 0.08 0. 10 TIME (SEC) P.7.5: S h o r t c i r c u i t s i m u l a t i o n . S i n u s o i d a l s o u r c e , z e r o a t t=0. New model v s . 750-Hz. F i r s t c y c l e s . 175 600 -i i • • • • • • • • • i 1 ' • • • ' • •> i — 1 — • <—<—>—•—'—<—<—i ' • — — > — i — < — < — i 0 0.010 0.020 0.030 0.040 TIME (SEC) P.7.6: Sho r t c i r c u i t s i m u l a t i o n . S i n u s o i d a l s o u r c e , peak a t t=0. New model v s . 60-Hz. 600 i 500 -j TIME (SEC) P.7.7: S h o r t c i r c u i t s i m u l a t i o n . S i n u s o i d a l s o u r c e , peak a t t=0. New model v s . 750-Hz. 176 6000 n 5000 o co Q L co 4000 5 3000 to ? 2000 I O O O H 1 1 I ' 1 0.0100 0.0200 TIME (SEC) 0.0300 0.0400 P.7.8: S h o r t c i r c u i t s i m u l a t i o n . Step e x c i t a t i o n . New model v s . 60-Hz. 3000 -> 2500 H o LO to 2000 S 1500 CO ? 1000 -| 500 A 0.0050 0.0100 TIME (SEC) i i 0.0150 • 0.0200 P.7.9: Sh o r t c i r c u i t s i m u l a t i o n . New model v s . 750-Hz. Step e x c i t a t i o n . P.7.10: Short circuit simulation. Step excitation. New model. Extended time response. -2.0 | i i i i 0 0.010 0.020 0.030 0.040 TIME (SEC) P.7.11: Open c i r c u i t s i m u l a t i o n . S i n u s o i d a l s o u r c e , peak at t=0. New model v s . 60-Hz. P.7.12: Open c i r c u i t s i m u l a t i o n . S i n u s o i d a l s o u r c e , peak a t t=0. New model v s . 750-Hz. 17 9 2.0 -, -0 .0 I '' i i i i 1 1 1 i i 0 0.005 0.010 0.015 0.020 0.025 0.030 TIME (SEC) P.7.13: Open c i r c u i t s i m u l a t i o n . Step e x c i t a t i o n . New model v s . 60-Hz. 2.0 -i -0 .0 I i i i i 0 0.005 0.010 0.015 0.020 0.025 0.030 TIME (SEC) P.7.14: Open c i r c u i t s i m u l a t i o n . Step e x c i t a t i o n . New model v s . 750-Hz. • 1 I I I ' I 0.005 0.010 0.015 0.020 TIME (SEC) P.7.15: Capacitive load. Step excitation. Load voltage. New model vs. 60-Hz. 2.5H TIME (SEC) P.7.16: Capacitive load. Step excitation. Load voltage. New model vs. 750-Hz. 300 4 0 - 1 — I — , — I — I — I — . — I — I — I — I — I — I -0.005 0.010 TIME (SEC) 0.015 P.7.17: Capacitive load. Step excitation. Load current. New model vs. 60-Hz. 300 n 250 -200 --100 --150 --200 --250 --300 I i i 0 0.005 0.010 0.015 TIME (SEC) P.7.18: Capacitive load. Step excitation. Load current. New model vs. 750-Hz. 182 0 . 5 i - 0 . 3 : - 0 . 4 : -0.5 1 i i -0 000 0.002 0.004 0.006 0.008 0,010 TIME (SEC) P.7.19: L-C l o a d . Step e x c i t a t i o n . Load v o l t a g e . New model v s . 60-Hz. 0.5 i 0.4 -- 0 . 3 : - 0 . 4 : -0 5 I • i • • i i i | | | [ i i i • i • i j -0 000 0 . 0 0 2 0 . 0 0 4 0 . 0 0 6 0 . 0 0 8 0 . 0 1 0 TIME (SEC) P.7.20: L-C l o a d . Step e x c i t a t i o n . New model v s . 750-Hz. Load v o l t a g e . 183 0.001 0.002 0.003 TIME (SEC) 0 .004 0.005 P.7.21: L-C load. Step excitation. Current in C. New model vs. 60-Hz. -0.000 0.001 0.002 0.003 TIME (SEC) 005 P.7.22: L-C load. Step excitation. New model vs. 750-Hz. Current in C. 184 TIME (SEC) P.7.23: L-C l o a d . Step e x c i t a t i o n . C u r r e n t i n L. New model v s . 60-Hz. 1600 -i -0.000 0.002 0.004 0.006 0.008 0.010 TIME (SEC) P.7.24: L-C l o a d . Step e x c i t a t i o n . C u r r e n t i n L. New model v s . 750-Hz. P.8.1(a): F i e l d test o s c i l l o g r a p h (BPA). S i n g l e - l i n e - t o -ground f a u l t on phase C. Receiving end voltages. P.8.1(b): Simulation of f i e l d test i n P.8.1(a) with the ne-l i n e model. (Only the zero sequence i s modelled as frequency dependent.) P.8.2: Comparison between v o l t a g e s a t phase b f o r : a) F i e l d t e s t o c i l l o g r a p h ; b) BPA's frequency-dependence s i m u l a t i o n ; c) New model s i m u l a t i o n . CONCLUSIONS The accurate modelling of transmission lines for time-domain electromagnetic transient calculations has been considered in this work. The transmission line model that has been developed is very simple, and yet, i t can very accurately represent the line behaviour over the entire frequency range of the signals. The main advantages of this new model as compared to previous frequency dependence formulations can be briefly summarized as follows: i) The accuracy of the model is not limited to certain specific frequencies, but includes the entire frequency range from the correct representation of dc levels to the highest frequency at which the system propagation function becomes negligible. It can, therefore, reliably model any possible system condition. i i ) The model is very efficient and does not appreciably increase the computer time of the transient simulations over the time required by constant-parameter representations. i i i ) The numerical routines for obtaining the parameters of the model are very general, numerically stable and efficient, and do not require user-defined control parameters. The overall process can then be completely automated without the need for user intervention other than to supply the physical configuration of the line. The generality of the principles involved in the development of the new line models presented in this work w i l l probably allow the extension of these principles to the modelling of frequency dependence and distributed losses in other system components, as for instance, power transformers, generators, and reactors. 187 Another area in which i t w i l l probably be possible to extend these concepts is in the modelling of multiphase systems with frequency-dependent transformation matrices. An important example of these systems is the case of underground cables. APPENDIX I SOME MATRIX RELATIONS^ I.1 Diagonalization of the Propagation Equations It is f i r s t assumed that the product (Z , Y , ) is diagonalizable by ph ph the matrix P whose columns are the eigenvectors of ) , that is P _ 1(Z ,Y ,)P = D (diagonal). (1.1) ph ph zy 6 Transposing both sides of (1.1), P C(Y I Z I ) (P "V = D 1 . ph ph zy Since Y , and Z . are symmetric and D is diagonal, then ph ph J zy Y t = Y ph ph Z t = Z ph ph D t = D zy zy and p t ( Y p h V ( P _ 1 ) t • Dzy • Let Q 1 = ? t, then eqn. 1.2 becomes Q.^Y ,Z >Q = D = D = D . (1.3) ph ph zy yz v ' Therefore, i t follows from eqns. 1.1 and 1.3 that i f the product (Z ,Y ,) ph ph is diagonalizable by a matrix P, (*) The mathematical demonstrations included in this Appendix have been provided by Luis Marti (M.A.Sc. candidate at U.B.C.) in consultation with Dr. L.M. Wedepohl (Dean, Faculty of Applied Science, U.B.C.). 189 190 P ' H z hY h)P = D (diagonal), (1.4) then, there exists a matrix Q such that the inverse product (Y , Z , ) is also ph ph diagonalizable, , Q" 1 ( Yph Z Ph ) Q = D ' ( I ' 5 ) where Of 1 = P t . (1.6) The eigenvector matrix P can be generalized as follows: P = P'D' , (1.7) where D' is an arbitrary, constant, diagonal matrix. This is possible because the eigenvectors of a matrix conserve their properties when multiplied by an arbitrary constant. Introducing equation 1.7 in equation 1.6, the following general relationship between the diagonalizing matrices P and Q is obtained: Cf 1 = (P'D')'' Q = (P'V^. , (1.8) where is again an arbitrary, constant, diagonal matrix (D = (D') . I.2 Diagonalization of the Drop Equations Given any two matrices A and B with non-zero diagonal elements and two diagonal matrices D^ and D^ with distinct elements, i f AB = D^ and BA = D^ , then D^ = a n& both A and B are diagonal. Proof: AB = D1 AD2 = ABA , 191 then AL>2 = D±A Let A - [a..] , » ± = [d k k] , etc. The product of any two matrices is given by [ a ] [ $ ] = [y] , where n y. . = E a., g, . . Then, i j k=l ik kj n n (AD0) . . = E a..d,.. = (D,A) . . = E d, .. a. . , 2 i j k=i lk 2kj 1 i j k=i l i k kj but since d, . = 0 for k ^ i , then kj (AD.).. = a..d,.. = (D.AJ.. = d, . . a.. , 2 i j i j 2 J J 1 l i ] I n 13 yielding a..d„ . . = a . .dn . . i j 2 J J i j I n a. .(d„. . - d, . .) = 0 i j 2 J J I n ' For i = j , a. .(d... - d. ..) = 0 , 11 2 n I n ' but since a.. 4 0 , then d„.. = d, .. and 11 2 n I n D l " D2 From this, and with i ^ j , a. .(d_. . - d„. .) = 0 , i j 2 J J 2n' which means that i f d... 4 d..., then a.. = 0 and therefore A is diagonal. 2 J J 211' i j 6 (Note that ^ ^ . j ^ ^ 2 i i * S a s u f f i c i e n t > D U t n o t necessary, condition.) If A 192 and are diagonal, i t follows directly from AB = that B is also diagonal. Q.E.D. Diagonality of the transformed matrices in the drop equations: Recalling equations 1.1 and 1.3, p ' 1 ( z p h V p = D Q~ 1 ( Yph zph ) Q = D • Introducing identity matrices on the left-hand sides, P" l zph (QCf 1 ) Yph P = D Q" l Yph ( P P" 1 ) Zph Q = D • and rearranging, Let ( P - ^ Q X Q - ^ P ) = D (1.9) (Q _ 1Y p hP)(P _ 1Z p hQ) = D . (I.10) A = p _ 1 V B = Q - ^ P , substituting in 1.9 and I.10, AB = D BA = D , and from the result proved f i r s t in this section, i t is concluded that i f D has distinct elements then A and B are diagonal; therefore, 193 P 1Z p hQ = D z (diagonal) (I.11) and Q _ 1Y P = D (diagonal) . (1.12) ph y In the particular case in which the elements of D are not distinct, that i s , the matrix (Z^Y^) ^ a s a t l e a s t t w o equal eigenvalues (e.g. positive and negative sequence modes in a perfectly balanced system), i t can be shown that an adequate selection of the eigenvector matrices P and Q can also (*) diagonalize Z ^ a n ^ ^pn» t n a t i s> satisfy eqns. I.11 and 1.12 (*) Dr. Wedepohl's lecture notes for course ELEC 552, Winter Session 1979-80, U.B.C., pp. 45-57. APPENDIX II INTEGRATION COEFFICIENTS II.1 Recursive Convolution Property Consider the following convolution integral; s(t) = /" f ( t - u ) e - a ( u - T ) d u . T At t-At, s(t-At) = f° f(t-At-u)e" a ( u" T )du Let v=At+u, then -(t-At) = £ f ( t _ v ) e - a ( v - T - A t ) d v ^ T+At or, rearranging and renaming the dummy variable v as u, s(t-At) = e a A t /" f ( t - u ) e - a ( u - T ) d u . T+At Separating the integral in equation II.1 into two parts, , , ..T+At.;,. . -a(u-T), , .<*> e . . -a(u-T), s(t) = / f(t-u)e du + I f(t-u)e du . X T+At Comparing the second integral with equation II.2 , . . -cxAt , , fT+At,, N -a(u-T), s(t) = e s(t-At) + / f(t-u)e 'du . T Evaluation of the one-time-step integral: The integral in equation II.3 can be evaluated as follows, Let v = u-T and = t-T, then 194 195 rT+Atc. , -a(u-T), r A t r / ^ - i - a v j / - r - r / \ / f(t-u)e du = / f ( t -v)e dv . (II.4) .j. o 1 • In a discrete process of solution, the function f is known only at discrete points: 0, At, 2At, .... To evaluate integral II.4, the intermediate values of f can be approximated by joining the known ends by a straight line, that i s , [f(t -At) - f ( t )] f(t x-v) = f ( t ; L ) + v , or f(t 1-v) = a + bv . Integral II.4 can now be evaluated: rAt .At, . -av, / = J (a+bv)e dv A . -av At . .At -av, re , a . n -aAt. J ae dv = a[ ] = —(1-e ) 0 -a J 0 a v ' At -av e~ a v A t / bve dv = b [ — — (-av-1) ] o a 2 o = \ [ l - e - a A t ( l + a A t ) ] Equation II.3 is fin a l l y given by s(t) = e- a A ts(t-At) + ( l - e - a A t ) ] f ( t - T ) a Ata z + I" + TS< l-e - O A t ) ] f (t-T-At) , (II.5) a A. to/ or s(t) = m s(t-At) + pf(t-T) + qf(t-T-At) , (II.6) (II.5) 196 where m, p, and q are constants. II.2 Dommel's Model for R-C Network Consider Dommel's discretized representation of an R-C combination [10] shown in Fig. II.1. e(t) •e(t) i ( t ) / r = At/2C Fig. II.1 From the equations in reference [10] , it,(t) is given by i h ( t ) = - i 2(t-At) - | | e(t-At) Solving the cir c u i t , but then i(t) = i ^ t ) + i 2 ( t ) , ±At) = , i 2 ( t ) = f e(t) + i h ( t ) = f e(t) - i 2(t-At) - f e(t-At) , i 2(t-At) = i(t-At) - i ^ t - A t ) = i(t-At) - -| e(t-At) , i 2 ( t ) = | | e(t) - i(t-At) + | e(t-At) - | | e(t-At) Combining terms and solving for e(t), 197 RA t In terms of the parameters of the partial fraction expansion of Z (s) eq (eqns. 4.13, 4.18, and 4.19), R = k/a and C = 1/k, i t is f i n a l l y obtained that . . rl-^(aAt)i , A . , . r At/2 -t • / \ e ( t ) - [l^Ktj] e ( t " A t ) + k [ l ^ ( a A t ) - ] l ( t ) + k [ l 4 ( a A t ) ] t ( t " A t ) > ("-a) or e(t) = m e(t-At) + p i(t) + q i(t-At) , (II.9) where m, p, and q are constants. BIBLIOGRAPHY [I] L.V. Bewley, Traveling Waves on Transmission Systems (book). Dover, New York, 1963 ( f i r s t published in 1933). [2] W.W. Lewis, The Protection of Transmission Systems Against Lightning (book). Dover, New York, 1965 ( f i r s t published in 1950). [3] E. Clarke, Circuit Analysis of A-C Power Systems (book, two volumes). Wiley, New York, 1950. [4] H.A. Peterson, Transients in Power Systems (book). Dover, New York, 1966 ( f i r s t published in 1951). [5] R. Rudenberg, Electrical Shock Waves in Power Systems (book). Harvard University Press, Cambridge, 1968 ( f i r s t published in German in 1962). [6] H.A. Peterson, "An Electric Circuit Transient Analyser". Gen. Elec. Rev., p. 394, 1939. [7] J.P. Bickford, N. Mullineux, and J.R. Reed, Computation of Power-System Transients (book). Peregrinus (for the IEE), Herts (England), 1976. [8] L.M. Wedepohl, "Application of Matrix Methods to the Solution of Travelling-Wave Phenomena in Polyphase Systems". Proc. IEE, Vol. 110 (12), pp. 2200-2212, December 1963. [9] D.E. Hedman, "Propagation on Overhead Transmission Lines: I. Theory of Modal Analysis; II. Earth-Conduction Effects and Practical Results", IEEE Trans., PAS-84, pp. 200-211, March 1965. [10] H.W. Dommel, "Digital Computer Solution of Electromagnetic Transients in Single-and Multiphase Networks". IEEE Trans., PAS-88, pp. 388-399, April 1969. [II] A. Budner, "Introduction of Frequency-Dependent Line Parameters into an Electromagnetic Transients Program". IEEE Trans., PAS-89, pp. 88-97, January 1970. [12] J.K. Snelson, "Propagation of Travelling Waves on Transmission Lines--Frequency Dependent Parameters". IEEE Trans., PAS-91, pp. 85-91, January/February 1972. [13] W.S. Meyer and H.W. Dommel, "Numerical Modelling of Frequency-Dependent Transmission-Line Parameters in an Electromagnetic Transients Program". IEEE Trans., PAS-93, pp. 1401-1409, September/October 1974. [14] A. Semlyen and A. Dabuleanu, "Fast and Accurate Switching Transient Calculations on Transmission Lines with Ground Return using Recursive Convolutions". IEEE Trans., PAS-94, pp. 561-571, March/April 1975. 1 9 8 199 [15] A. Ametani, "A Highly Efficient Method for Calculating Transmission Line Transients". IEEE Trans., PAS-95, pp. 1545-1551, September/ October 1976. [16] A. Semlyen and R.A. Roth, "Calculation of Exponential Step Responses -Accurately for Three Base Frequencies". IEEE Trans., PAS-96, pp. 667-672, March/April 1977. [17] A. Morched and A. Semlyen, "Transmission Line Step Response Calculation by Least Square Frequency Domain Fitting". IEEE PES Summer Meeting, Portland, OR, July 1976. [18] L.F. Woodruff, Principles of Electric Power Transmission (book). Wiley, New York, 1938. [19] J.R. Carson, "Wave Propagation in Overhead Wires with Ground Return". Bell Syst. Tech. J., vol. 5, pp. 539-554, 1926. [20] F. Pollaczek, "Uber das Feld einer unendlich langen wechselstromdurch-flossenen Einfachleitung"("About the f i e l d of an i n f i n i t e ac excited single-phase line".) Elek. Nachrichtentech, vol. 3, pp. 339-359, 1926. Translation by J.B. Pomey into French, "Sur le champ produit par un conducteur simple infiniment long parcouru par un currant alternativ". Gen. Elec. Rev., vol. 29, pp. 851-867, 1931. [21] H. Karrenbauer, "Ausbreitung von Wanderwellen bei verschiedenen Anordnungen von Freileitungen im Hinblick auf die Form der Einschwing-spannung bei Abstandskurzschlussen". ("Propagation of travelling waves for various overhead line tower configurations, with respect to the shape of the transient recovery voltage for short line faults".) Ph.D. dissertation, Technical University Munich, Germany, 1967. [22] P.C. Magnusson, "Travelling Waves on Multi-Conductor Open-Wire Lines — A Numerical Survey of the Effects of Frequency Dependence of Modal Composition". IEEE Trans., PAS-92, pp. 999-1008, May/June 1973. [23] R.G. Wasley and S. Selvarinayagamoorthy, "Approximate Frequency-Response Values for Transmission-Line Transient Analysis". IEE Proc, vol. 131 (4), pp. 281-286, April 1974. [24] H.W. Dommel and W.S. Meyer, "Computation of Electromagnetic Transients". IEEE P r o c , vol. 62(7), pp. 983-993, July 1974. [25] E. Groschupf, "Simulation transienter Vorgange auf Leitungssystemen der Hochspannungs-Gleichstrom-und-Drehstrom-Ubertragung". ("Simulation of transient phenomena in systems with HVDC and ac transmission lines".) Ph.D. dissertation, Technical University Darmstadt, Germany, 1976. [26] UBC - Overhead Line Parameters Program. Originally written by H. W. Dommel at the Bonneville Power Administration, Portland, Oregon. Modified at the University of British Columbia (UBC), Canada, by I. I. Dommer,. K.C. Lee., and T. Hung. User's manual, UBC, August 1980. 200 [27] H.W. Dommel, Notes on Power Systems Analysis. University of British Columbia, Department of Electrical Engineering, 1975. [28] E.A. Guillemin, Synthesis of Passive Networks (book). Wiley, New York. 1957. [29] S. Karni, Network Theory: Analysis and Synthesis (book). Allyn and Bacon, Boston, 1966. [30] A. Semlyen, "Contributions to the Theory of Calculation of Electro-magnetic Transients on Transmission Lines with Frequency Dependent Parameters". IEEE Trans., PAS-100, pp. 848-856, February 1981. [31] S.W. Director, "Survey of Circuit-Oriented Optimization Techniques". IEEE Trans., CT-18, pp. 3-10, January 1971. [32] G. Szentirmai, "Computer Aids in F i l t e r Design: A Review". IEEE Trans., CT-18, pp. 35-40, January 1971. [33] H.W. Bode, Network Analysis and Feedback Amplifier Design (book). Van Nostrand, New York, 1945. [34] A. Papoulis, The Fourier Integral and i t s Applications (book). McGraw-Hill, New York, 1962. [35] B.W. Arden (editor), What Can Be Automated? The Computer Science and Engineering Research Study (COSERS) (book). The MIT Press, Cambridge, Massachusetts, 1980. [36] R.W. Hornbeck, Numerical Methods (book). Quantum, New York, 1975. [37] H.W. Dommel and N. Sato, "Fast Transient Stability Solutions". IEEE Trans., PAS-72, pp. 1643-1650, July/August 1972. [38] UBC - Electromagnetic Transients Program (EMTP). Originally written by H.W. Dommel at the Bonneville Power Administration, Portland, Oregon. Modified at the University of British Columbia (UBC), Canada, by H.W. Dommel. User's manual, UBC, August 1978 ( f i r s t published in 1976) . [39] J. Umoto and T. Hara, "A New Digital Analysis of Surge Performances in Electric Power Networks Utilizing the Convolution Integral". (Trans-lated from Denki Gakkai Zasshi, vol. 91(5), pp. 907-916, May 1971.) Electrical Engineering in Japan, vol. 91(3), pp. 48-57, 1971.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The problem of frequency dependence in transmission...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The problem of frequency dependence in transmission line modelling Martí, José R. 1980
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | The problem of frequency dependence in transmission line modelling |
Creator |
Martí, José R. |
Publisher | University of British Columbia |
Date Issued | 1980 |
Description | In this work, the accurate representation of transmission lines for the digital simulation of electromagnetic transients in power systems has been examined. A model has been developed that accounts for the frequency dependence and distributed nature of the line parameters over the entire frequency range. This model can easily be incorporated into a time-domain network solution of the complete power system. The model consists simply of a constant resistence in parallel with a current source evaluated at each time step of the solution. The equivalent resistance results from a finite-step-width discretization of the differential equations of a resistance-capacitance (R-C) network that simulates the line characteristic impedance. The equivalent current source accounts for the time delays and attenuations of the different frequency components of the travelling waves and for the discretization of the time-domain equations. Rational-function approximations are used to synthesize the R-C network and the line propagation ("weighting") function in the frequency domain. These rational approximations allow the corresponding time-domain functions to be obtained directly in a closed-form, thus circumventing the need for numerical inverse Fourier transformations. The numerical technique used to obtain the rational functions yields very accurate, high-order approximations. This technique is based on a direct, step-by-step allocation (and reallocation) of poles and zeros and avoids the instability problems which can be encountered with optimization techniques based on search methods. A series of analytical evaluations and simulation tests were performed in order to assess the validity of the model. The results of these tests show that the model is accurate, fast, and reliable. The model was incorporated into the code of the University of British Columbia's version of Dr. H.W. Dommel's Electromagnetic Transients Program (EMTP). i |
Subject |
Transients (Electricity) Electric lines --Models |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-03-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0095571 |
URI | http://hdl.handle.net/2429/23154 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1981_A1 M37_3.pdf [ 7.98MB ]
- Metadata
- JSON: 831-1.0095571.json
- JSON-LD: 831-1.0095571-ld.json
- RDF/XML (Pretty): 831-1.0095571-rdf.xml
- RDF/JSON: 831-1.0095571-rdf.json
- Turtle: 831-1.0095571-turtle.txt
- N-Triples: 831-1.0095571-rdf-ntriples.txt
- Original Record: 831-1.0095571-source.json
- Full Text
- 831-1.0095571-fulltext.txt
- Citation
- 831-1.0095571.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0095571/manifest