SIMULATION O F LIGHTNING SURGES O N TRANSMISSION LINES By HUYEN V. NGUYEN B.Sc, Elec. Engr., Gannon University, 1987 M.Engr., Electric Power Engr., Rensselaer Polytechnic Institute, 1988 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES D E P A R T M E N T OF E L E C T R I C A L ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C O L U M B I A February 1996 © Huyen V. Nguyen, 1996 In presenting degree this thesis in at the University of partial fulfilment British Columbia, of the requirements for an I agree that the Library shall make it freely available for reference and study. I further agree that permission copying of department this or thesis by his for scholarly or her purposes may representatives. be granted It is for extensive by the head understood that publication of this thesis for financial gain shall not be allowed without permission. Department of ZieCTXlCAL The University of British Columbia Vancouver, Canada Date DE-6 (2/88) h^L re, wc, advanced €M£IlV€£fZlh)fr- of my copying or my written ii ABSTRACT Computer simulations with programs such as the EMTP are widely accepted for lightning transient studies, to design the system in such a way that equipment failures and transmission line outages are prevented. In these studies, an accurate representation of the transmission lines and towers is required. Two new and improved time-domain models for the representation of these components are presented in this thesis. They are the frequency-dependent and the exponential transmission line models. The former is based on synthesizing the line functions directly in phase co-ordinates, and is used to represent multiphase transmission lines, while the latter is a single-phase model, which also synthesizes the line functions, to be used for tower representations. These models are valid for the required frequency range of 1 Hz to 1 MHz. They can be used for lightning studies as well as for many other types of transient studies. Extensive time-domain simulations have been performed to compare the results with those produced by existing line models. The results indicate that the new models are more accurate than some of the existing models, and numerically stable. Simulations of typical lightning surges have also been performed to compare the solutions produced by the developed models with actual measurements. With these new models, simulations match field recorded transient waveforms very well. In addition, a simple, fast and accurate method based on the Newton-Raphson iteration technique is introduced for finding transmission line modal parameters. It overcomes the eigenvector and eigenvalue switching problems experienced by most other methods. With a proper constraint equation for scaling the eigenvectors, this proposed method produces modal parameters which can reasonably be fitted with rational functions of minimum-phase type. With this, the modal transmission line modelling can be improved. The routine is very stable and can be used to diagonalize the Y Z product for uncommon lines with strong asymmetry. It has also been used to obtain the phase-domain parameters for the new frequency-dependent line model. iii Table of contents Abstract ii Table of contents List of Tables iii — . : •. — v List of Figures v Acknowledgement ix C H A P T E R 1: Introduction 1 1.1) 1.2) 1.3) 1.4) 1.5) Background Types of lightning transients on overhead lines Lightning impacts on high voltage power systems Thesis outline Thesis contributions 1 2 2 4 5 C H A P T E R 2: Literature Review 2.1) 2.2) 2.3) 2.4) 12 Lightning waveform measurements in power systems Lightning stroke representation in engineering calculations Tower surge impedance and its response to lightning strokes Transmission line characteristics for lightning studies C H A P T E R 3: Modelling of Single-Phase Nonuniform Transmission Lines for Representation of Transmission Towers 3.1) Overview 3.2) Exponential transmission line 3.3) Synthesis of exponential line functions 3.4) Reflection factors 3.5) E M T P simulations : — C H A P T E R 4: Frequency-Dependent Transformation Matrices for Untransposed Transmission Lines using Newton-Raphson Method 4.1) 4.2) 4.3) 4.4) Overview Equations for multiphase transmission lines — Eigenproblem solution with Newton-Raphson (NR) method Modal parameters for multiphase transmission lines from Newton-Raphson method 4.5) Open/Short -circuit tests in the frequency domain 4.6) Validation of the Newton-Raphson algorithm with the power method 12 13 15 17 20 20 21 25 31 35 37 37 38 ...42 45 : — 54 60 C H A P T E R 5: Frequency-Dependent Overhead Transmission Line Model in the Phase Domain 5.1) 5.2) 5.3) 5.4) 5.5) 5.6) Overview Phase-domain transmission line modelling Synthesis of the propagation and characteristic matrices Synthesis of Y and A for asymmetrical transmission line configurations Open/Short -circuit tests in the frequency domain Time-domain simulation tests c C H A P T E R 6: Simulation of Lightning Surges and Field Test Comparisons 6.1) 6.2) 6.3) 6.4) 6.5) Overview Tower/ground wire surges Simulation of a recorded shielding failure event in a 765 kV system Shielding-failure flashover Back flashover 62 62 62 65 73 77 87 91 91 91 96 103 104 C H A P T E R 7: Future Research in Transmission Line and Transmission Tower Modelling Ill C H A P T E R 8: Conclusions 113 APPENDIX A: Asymptotic Behaviour of Modal Transformation Matrices at Very Low and High Frequencies 115 APPENDIX B: Recursive Convolutions ..118 APPENDIX C: Interfacing User-defined Models with the EMTP version MicroTran . .121 APPENDIX D: Estimation of Maximum Lightning Shielding Failure Stroke on a Typical 765 kV Line Based on Electrogeometric Theory 124 APPENDIX E: Representation of a Long Line with Linear A C Steady-State Initial Conditions 128 REFERENCES 130 V List of Tables Table 5.1: Average number of poles allocated for synthesized functions 74 Table 5.2: Steady-state induced voltages (RMS) 90 List of Figures Figure 1.1- Multiple shielding failures waveforms recorded on March 26, 1991 6 Figure 1.2 - Expanded view of stroke #1 of Figure 1.1 7 Figure 1.3 - Expanded view of stroke #2 of Figure 1.1 8 Figure 1.4 -Expanded view of stroke #3 of Figure 1.1 9 Figure 1.5 - Shielding failure waveforms recorded on April 12, 1991 10 Figure 1.6- Expanded view of strike location of Figure 1.5 11 Figure 3.1 - Single-phase exponential line 21 Figure 3.2 - Equivalent circuit for the exponential line 25 Figure 3.3 - Magnitude of Z 27 Figure 3.4 - Phase of Z c k in Q c k in degrees — 27 Figure 3.5 - Magnitude of A 28 Figure 3.6 - Phase of A in radians 28 Figure 3.7 - Propagation function a(t) 29 Figure 3.8 - Magnitude of Z ' 30 Figure 3.9 - Phase of Z ' cm cm in Q in degrees 30 Figure 3.10- Single-phase exponential line terminated at node m with Z L 31 Figure 3.11 - Magnitude of Z in Q from Eq. (3.24) 33 Figure 3.12- Magnitude of Reflection factor for Z = 150 Q 34 Figure 3.13- Magnitude of Reflection factor for Z = 100 Q 34 Figure 3.14 - Receiving end open-circuit voltage (p.u.) 36 Figure 3.15- Receiving end short-circuit current (p.u.) 36 Figure 4.1 - Multiphase transmission line system 39 Figure 4.2 - Vertical double-circuit line 47 in L L VI Figure 4.3 - Elements of column 1 of the matrix T , 50 Figure 4.4 - Elements of the matrix Y C M 51 Figure 4.5 - Elements of the matrix A M 52 Figure 4.6 - Phase angle of element A . , calculated from Eq. (4.21) 53 Figure 4.7 - Phase angle of element A _, calculated with atan function 53 Figure 4.8 - Receiving end open-circuit voltages (p.u. magnitude) 56 Figure 4.9 - Receiving end short-circuit currents (p.u. magnitude) 57 Figure 4.10 - Receiving end open-circuit voltages (p.u. magnitude) 58 Figure 4.11 - Receiving end short-circuit currents (p.u. magnitude) 59 Figure 4.12- Locations of largest elements of the correlation matrix 61 Figure 5.1 - Multiphase transmission line equivalent circuit in the time domain 65 M M Figure 5.2 - Elements of column 1 of matrix Y for a double-circuit line 69 Figure 5.3 - Elements of column 1 of matrix P 70 c Figure 5.4a - Phase angle (degrees) of elements P and P before multiplying with -1.0 . 71 21 Figure 5.4b - Phase angle (degrees) of elements A 3I and A n 2 1 from Eq. (5.12) 71 Figure 5.5 - Time-domain propagation functions: a (t), -a (t), and -a (t) 72 Figure 5.6 - Current source with matching admittance Y (co) connected to node m 73 Figure 5.7 - Railroad near a three-phase 138 kV line 74 u 21 31 c Figure 5.8 - Elements of column 1 of matrix Y for an asymmetrical line c Figure 5.9 - Elements of column 1 of matrix A for an asymmetrical line 75 76 Figure 5.10 - Receiving end open-circuit voltages (p.u. magnitude) for line of Figure 4.2 79 Figure 5.11 - Receiving end short-circuit currents (p.u. magnitude) for line of Figure 4.2 80 Figure 5.12 - Receiving end open-circuit voltages (p.u. magnitude) for line of Figure 5.7 81 Figure 5.13 - Receiving end short-circuit currents (p.u. magnitude) for line of Figure 5.7 82 Figure 5.14- Receiving end open-circuit voltages (p.u. magnitude) for line of Figure 4.2 83 Figure 5.15 - Receiving end short-circuit currents (p.u. magnitude) for line of Figure 4.2 84 Figure 5.16 - Receiving end open-circuit voltages (p.u. magnitude) for line of Figure 5.7 85 Figure 5.17 - Receiving end short-circuit currents (p.u. magnitude) for line of Figure 5.7 86 Figure 5.18 - Receiving end phase 1 open-circuit voltage (p.u.) 87 Vll Figure 5.19 - Magnetically induced voltages on rails 89 Figure 5.20 - Magnetically induced voltage on rail #1 89 Figure 6.1 - Circuit for EMTP simulations 93 Figure 6.2 - Conductor geometry at average height 94 Figure 6.3 - Tower representation -94 Figure 6.4 - Simulation of Ishii et al experiment and measurement from [18] 95 Figure 6.5 - Simplified diagram of 765 kV system surrounding Rockport 99 Figure 6.6 - Lightning current source (A) 99 Figure 6.7 - Lightning shielding failure waveforms recorded on Oct. 7, 1990 100 Figure 6.8 - Expanded view of recorded shielding failure waveforms 101 Figure 6.9 - E M T P simulations of the shielding failure event in Figures 6.7 & 6.8 102 Figure 6.10 - E M T P simulations of a shielding-failure flashover event 107 Figure 6.11 - E M T P circuit for simulation of a back flashover event 108 Figure 6.12 - Lightning current source (A) 108 Figure 6.13 - Simulated voltages at tower #7 with a 15 Q footing resistance 109 Figure 6.14 - Simulated voltages at tower #7 with a 180 Q footing resistance 110 Figure C l - Equivalent circuit seen by the transmission line terminals 121 Figure D l - Shielding angle of the right conductor bundle 124 Figure D2 - Incomplete Shielding, Width X is uncovered 127 Figure D3 - Effective shielding, Unprotected width X is reduced to zero 127 Figure E l - Single transmission line 128 Figure E2 - Equivalent circuit of a long line without initial conditions 129 Figure E3 - Equivalent circuit of a long line with linear A C initial conditions 129 s s viii In memories of Reverend Monsignor Wilfrid J. Nash and my two dear sisters Dung and Xuan ix Acknowledgement I would like to express my deepest appreciation to all the people who contributed directly or indirectly to this thesis project. The following thanks are expressed specifically to those who have and had enabled me in realizing my full potential: To my thesis supervisors, Drs. Hermann Dommel and Jose Marti, whom I respect and admire tremendously. Not enough can be said about my appreciation for all they have helped in giving guidance, encouragement, as well as arranging much-needed financial support. Their unique, patient, and reasonable ways of getting ideas across had allowed me to learn and perform so freely. To Dr. Martin Wedepohl for his patience and effort in introducing the Newton-Raphson method of Chapter 4. To Dr. Luis Marti for his help and support as Industrial Technical Advisor. To Drs. A.J.F Keri, J.M. Schneider and C.H. Shih for providing field recordings and system parameters for the shielding failure and railroad cases, and also for training and guiding me in conducting assigned projects during my stay at the American Electric Power Service Corporation, Columbus, Ohio. It was during this period that I have become so interested in what is called "power system transients" To Dr. Jerry Selvaggi for passing his love of engineering on to me. To all my colleagues, especially to Fernando Castellanos and Luis Linares, in the power group at U B C for creating such a warm and friendly research environment. The financial assistance of the Natural Sciences and Engineering Research Council of Canada through a grant to Dr. Hermann Dommel, and the graduate fellowships awarded by the University of British Columbia are gratefully acknowledged. I am also deeply grateful to my dearest uncle Van for his encouragement and help in liberating me to America. Thanks are addressed to all my brothers Doan, Tuan , Huy, Van, and to my oldest sister Hang, for their help and support. Thanks are also given to all my friends at the Vietnamese Catholic Church choir (1992-1996) in Vancouver, B.C., for their kindness in many ways. It is through this organization that I have been blessed to have known my wonderful wife Hien. Thanks are deeply expressed for her love, patience and support during these years. And finally, I would like to dedicate all my achievements to my dear parents. Thanks for their love and support through the years. 1 CHAPTER 1 Introduction •i , . . • . 1.1) Background Analyses of field recording events over many years [2,8,9] have suggested that high overvoltage transients caused by lightning are considered to be the chief source of disturbances in high voltage transmission systems. Lightning phenomena have been studied throughout this century, although the process that initiates lightning is still not completely understood. This has led to the development of various theories in attempting to describe the so-called lightning-mechanism. There is, however, a consensus that lightning starts from the charge separation process (positive and negative), which is due to the transportation of lightweight particles to higher regions by the rapid updrafts of moist air, usually in hot humid areas. This vertical charge separation is known as the vertical thunderstorm dipole. It can be formed within the cloud or between the cloud and the earth, which creates electric fields that eventually bring about the electric breakdown known as lightning [3,4]. It is not within the scope of this thesis to fully explain the details of how the lightning phenomenon develops. What concerns the power system engineer is the severity of overvoltages caused by lightning in the transmission grid. The overvoltages introduced by lightning have traditionally been estimated using conventional and simplified methods [3,5,7]. More involved calculations have become possible with digital computer programs such as the 2 Electromagnetic Transients Program (EMTP) [1]. In such programs, each power system component can be modelled in great detail. 1.2) Types of lightning transients on overhead lines The characteristics of lightning surges on overhead transmission lines, which result from lightning strokes, depend on how they are caused. They can be broadly divided into four types (ignoring the induction case): • a) Tower/ground wire surges - The stroke terminates on the tower structure/ground wires without any back flashover to the phase conductors. • b) Shielding failure - The stroke passes through the protective zone of the ground wires and terminates on the phase conductor(s). • c) Back flashover - The same as (a), but followed by a flashover to the phase conduc- tors). • This type of flashover is called back flashover. d) Shielding-failure flashover - The same as (b), but followed by a forward flashover to the ground/ground wires or tower. 1.3) Lightning impacts on high voltage power systems Lightning overvoltages in power systems are either caused by a lightning stroke to the line conductors, the tower structure and the ground wires, or they can also be created by an induction process from a nearby stroke to ground. The latter causes insignificant overvoltages in high voltage systems, but the former type may result in system faults from either forward or back flashovers. These flashovers can damage terminal equipment. The 3 effects produced by those events which do not cause flashovers can also be severe. There have been many interesting shielding-failure events captured in the AEP (American Electric Power) 765 kV system, which created moderate overvoltages [2]. Figure 1.1 depicts the phase-to-ground waveforms at a 765 kV transformer terminal, which resulted from a multiple shielding-failure event, captured during heavy thunderstorms. This event is very unique since three consecutive different strokes, within a 110-ms period, terminated directly on the phase A conductor approximately 149 miles away; flashover occurred on the last stroke. The expanded views at the stroke locations are exhibited in Figures 1.2-1.4. The time between the first and the second, and the second and third strokes are approximately 65 and 45 ms, respectively. Figure 1.5 also shows another extraordinary shielding-failure event, which occurred near the negative peak of the system 60 Hz voltage. The expanded view is shown in Figure 1.6. Phase C is the struck phase in this case; it registered a maximum magnitude of 1.80 p.u.! These overvoltages, with non-standard waveshapes, are fully seen by the system equipment at the substations. The resulting effects of these non-standard waveforms impinging across the substation transformers, for example, are not yet well known to the transformer manufacturers. Many unexplained E H V transformer failures [2] may partially have been caused by these effects. Thus, it is important to realistically simulate and fully understand all the possible events created by lightning in power systems. The ensuing outcome can then be helpful for the future development of improved protection schemes. 4 1.4 Thesis outline Power system transients, especially those resulting from lightning, are by and large computed with general-purpose programs such as the EMTP [44]. The goal of this thesis is to develop new and improved transmission line and transmission tower models for lightning surge simulations for use in such computer programs. Overhead lines and transmission towers are both entirely exposed to lightning activities, and therefore their behaviour influences greatly the travelling waveforms reaching the substations. For modelling other system components, guidelines are available which produce results with reasonable accuracy [6], [7], and [44]. The thesis gives a brief literature review in Chapter 2, where the lightning return-stroke current measurements and models for overhead transmission lines and towers are discussed. A nonuniform transmission line model, based on an exponential variation of the line parameters, is presented in Chapter 3. This model is intended for the representation of transmission towers. For developing a time-domain frequency-dependent transmission line model, an eigenvector/eigenvalue routine based on the Newton-Raphson technique is investigated and validated in Chapter 4. A phase-domain frequency-dependent transmission line model is described in Chapter 5, along with various frequency and time domain test results. In Chapter 6, simulation results of common lightning events, along with field test comparisons, are presented. Possible for future research are outlined in Chapter 7, and conclusions about the research findings of the thesis are made in Chapter 8. 5 1.4 Thesis contributions The original contributions from the work in this thesis are: • Implementing and validating the Newton-Raphson method for evaluating the modal parameters of overhead transmission lines as smooth functions of frequency [37]. • A time-domain algorithm for modelling single-phase nonuniform transmission lines for the representation of transmission towers [22, 46]. • A time-domain algorithm for direct phase-domain modelling of frequency-dependent overhead transmission lines [47]. • Modelling techniques for simulation of typical lightning events. 6 7 at o > o Q. O O a. o at o > m a> in in t o Q. o o a: O) 2 o > - 6 to a O o CC "3~n3 3 0.5 ir^r Time (ms) 3 1.5 Figure 1 . 2 - Expanded view of stroke # 1 of Figure 1 . 1 x - 2 . 52. - 3 1 0 5 Time (ms) . > - 3 . 5 . Time (ms) Figure 1.4 -Expanded view of stroke #3 of Figure 1.1 10 12 CHAPTER 2 Literature Review 2.1) Lightning waveform measurements in power systems The overvoltages caused by lightning in power system are influenced by various characteristics of the lightning flashes. These characteristics include the stroke leader mechanism, the length-of-path, current polarity, time duration, charge distribution in the stroke channel and the total number of strokes. Many researchers have attempted to measure lightning currents from strokes to high towers. For example, the C N tower in Toronto is equipped with instrumentation to measure lightning stroke currents. Others have used lightning-flash counters to obtain statistics of the frequency of lightning strokes to earth, while others have used the power system itself as a means to record lightning transients, which impinge upon terminal equipment. Recordings from a monitoring system, which has been installed in the 765 kV system of American Electric Power Service Corporation, will be used in this thesis to compare simulation results against field measurements. Each type of transient phenomenon occurring in a power system requires a specific measuring technique, which depends primarily on the frequencies involved in the phenomenon. Lightning phenomena are of the fast transient type (10 kHz - 3 M H z [6]). Its waveforms in power systems have been captured with various monitoring systems [2,8,9,10]. Such monitoring systems must have a wide frequency bandwidth and fast sampling rates. 13 The system described in [2] captures lightning events in the form of voltages, through voltage dividers connected to transformer bushing capacitance taps. This monitoring system is of particular interest to the author since he managed and maintained it for a three-year period during his stay at the American Electric Power Service Corporation. During this time, it captured many interesting lightning events, such as those shown in Figures 1.1-1.6. 2.2) Lightning stroke representation in engineering calculations One of the most critical factors in lightning surge simulations is the representation of the • lightning flash itself. It has been accepted that when lightning strikes a power line, the lightning stroke injects a current source into the power system [3,7]. The injected current representing the lightning return-stroke is generally modelled as a current source, with the stroke-channel impedance often being ignored because it is relatively large. A single lightning flash may contain several strokes, with the average being three. If each stroke has a certain waveshape and magnitude, it may be necessary to represent the flash as a series of current sources in simulations. To understand the significant effects created by a multiple shielding-failure event, this kind of multi-stroke representation with lightning current sources is needed. In most cases, however, a single current source for the representation of only the leader stroke is justified [7]. Some other factors may be important as well. For example, Greenwood [3] mentions that a more accurate representation of lightning must take the leader stroke and the prestrike into account. Some researchers have theorizing that the effect of charge distribution in the stroke-channel should be included in lightning calculations too. Researchers have held different viewpoints about this particular parameter. 14 Wagner-Hileman [11] stated that the effect of the charge distribution in the stroke-channel on the line insulator voltage is as important as the effect of the injected current. Choy and Darveniza [12] re-examined this theory and confirmed the significant influence of this parameter on the predicted flashover rates. In their analysis, they assumed that the length of an upward streamer is one half of the striking distance. They found that the resulting insulator voltage due to charge in the stroke-channel decreases significantly as the length of the streamer increases. Subsequently, Darveniza et al [12] reversed this opinion and recommended to ignore this 'presumably' important parameter in backflash calculations. Most researchers today seem to model only the effects of the injected current stroke. This approach is also adopted in this thesis. Thus, the effect of charge in the stroke-channel will not be considered here. As far as the representation of the injected lightning current is concerned, factors such as the waveshape and the rate-of-rise significantly influence the overvoltages which lightning can create. Most earlier studies assumed the return-stroke current as a triangular-ramp function, because of the need for simplicity in the analyses involved. IEEE recommendations for surge arrester studies still use such simple waveforms. With the advent of digital simula- tions, which make more complicated source representations easy, more recent studies often used a double exponential current source. The most recent lightning waveform, proposed by O G R E Working Group 33 [7], gives a somewhat more realistic picture, since it has the concave wavefront, similar to that observed in measurements. It is therefore recommended that this new waveform be used in digital simulations. The latest U B C E M T P version Micro-. Tran has this function implemented. 15 2.3) Tower surge impedance and its response to lightning strokes The overvoltage created by a lightning stroke depends on the surge response characteristics of the object it strikes. Transmission towers are difficult to represent in lightning surge simulations. From a wave propagation standpoint, the tower is a non-uniform structure, although many models in the literature assume that it is uniform. Some non-uniform models have been proposed [19-24], but the mathematics involved in most of these models are complicated and not easy to implement in general purpose programs such as the Electromagnetic Transients Program (EMTP). With the EMTP, transmission towers have been modelled as pure resistances, as pure inductances, or as lossless constant-parameter transmission lines, but the accuracy of such representations has seldom been verified. This has motivated recent research on more accurate tower models [18]. It is relevant to mention that poor lightning performance has been reported on some transmission lines, even if they had very good shielding and low footing resistances [27]. Some lines with low footing resistances have actually displayed outage rates as high as the ones with high values [28]. These unexpected results may well be caused by wave propagation effects up and down the tower. For example, Chisholm [28] pointed out that the influence of the transmission tower on lightning performance of transmission lines may well have been overlooked in many lightning surge simulations. In the past, it has often been said that it is not practical to model transmission towers in detail, because there are already so many uncertainties involved in lightning surge simulations (see, for example, the discussions to [17]). Recent line performance data indicate that the actual lightning outage rates may be underestimated with simplistic models. It can also 16 be shown that the wave propagations up and down the tower after a back flashover have an important influence on the waveform which reaches the equipment in the substation. Better tower models are therefore needed. The transmission tower is a non-uniform transmission line structure, whose surge impedance (Z) varies as the" surge travels through it. Wagner and Hileman [13] initially derived and proposed a cylindrical tower model with the equation, Z=601n 2lf (2.1) where t is the travel time counted from the instant the current has entered the tower top, c is the velocity of light and r is the radius of the cylinder. Equation (2.1) illustrates that the tower surge impedance is not constant; it is lowest at the top and increases as the wave progresses down to the base. Kawai [14] performed measurements on isolated towers (without ground wires connected) and obtained a trend similar to the Wagner-Hileman equation, although the magnitudes were appreciably lower. However, the change in surge impedance in [13] and [14] is opposite to that found in most subsequent experiments [15]. Recognizing the inconvenience of the time-variation in equation (1), Sargent and Darveniza [16] developed a conical tower formula with a constant surge impedance, Z=601n (2.2) where S is the sine of the half-angle of the cone. Chisholm et al [17] later performed experiments as well as analytical calculations, and discovered that the position-dependent characteristics of the tower surge impedance are 17 dependent upon the direction of the injected current source. For a vertical stroke (current terminating at the tower top), the trend follows that of Kawai's and Wagner-Hileman's results. The trend is the opposite for the horizontal situation - where the stroke current terminates somewhere on the ground wire and progresses horizontally as a travelling wave towards the tower. To accommodate the case of the horizontal wave, they devised another equation which characterizes that of an inverted cone, (2.3) where 0 is the half-angle of the cone. Ishii et al [18] recently performed excellent low current measurements and ascertained that the presence of the connected ground wires alters the change of the measured impedance. Their results confirmed that the surge impedance is higher at the top and lower at the bottom when the tested tower has ground wires connected to it. These results will be used to verify the proposed nonuniform tower model, which is described in detail in Chapter 3. 2.4) Transmission line characteristics for lightning studies The modelling of transmission lines in digital simulations is well described in the literature; e;g., [30] to [41], to mention just a few. Most studies in the past have relied heavily on single-phase representations.: Due to. major improvements of the line models in the last two decades, detailed.multi-phase line models, are now used more and more. However, despite these improvements, limitations still exist. There are some guidelines for choosing the existing line models for typical Ughtning studies [5, 7, 31, 32, ...], but further improvements are 18 still required in order to take into account the complete frequency-dependent nature of untransposed overhead transmission lines. There are presently four main overhead transmission line models available in the EMTP, namely, the nominal Pi-circuit, the exact Pi-circuit^ the constant distributed parameter, and the frequency-dependent distributed parameter models [31], [36, pp. 4-66 to 4-107]. The assumptions used in these line models are clearly described in the EMTP Theory Book [36]. Both Pi-circuit models cannot be used for lightning or fast transient studies because they are valid only for steady-state or low-frequency transient analysis. If the constant-parameter model is used in lightning studies, the line should usually be modelled as untransposed, and the transformation matrix for transforming phase to mode quantities should be evaluated in the hundreds of kHz range. Unfortunately, the option of lumping R in a few places, which is used in many studies to take the line losses into account, is not usable for long lines in lightning studies, since the total line series resistance at high frequencies is no longer much less than the line surge impedance. This is one of the important restrictions which the con- stant-parameter line model has. Ignoring the line losses altogether might be a legitimate approximation for lightning studies, but this implies that all modes travel at the same speed of light. In reality however, the ground mode does travel somewhat slower (-5%) than the aerial modes in the lightning frequency range. If the studies are concerned with short line distances (a few tower spans), this model is considered to be accurate enough, since the mode separation over such a short distance is not noticeable [36, p. 4-57]. A better high frequency lossless constant-parameter model has been suggested [42], which starts from a line parameter calculation (possibly at 500 kHz) with losses included, but the series resistances are 19 set to zero just before the modal parameters are calculated. This gives a realistic wave velocity for the ground mode. One limitation should be kept in mind, though, namely that the resistance for the ground return mode in lightning transients is highly frequency dependent [32]. For this reason, it is best to use the frequency-dependent model. However, even this model has the limitation of using constant transformation matrices. Constant transformation matrices are usually accurate enough for single and flat circuit lines. They are less accurate for double-circuit lines and probably worse for triple circuit lines. In general, the constant transformation matrix assumption is questionable if the asymmetry of the line conductor configuration is strong [33]. Thus, a full wideband frequency-dependent transmission line model is desired [43]. This line model should provide correct answers in both steady-state and transient studies. This main objective has been achieved in this thesis. The modelling details of this new transmission line model are presented in Chapter 5. It is worth mentioning that a full frequency-dependent modal-domain model was developed earlier to achieve the same goal [35], but primarily for underground cables. In parallel with the work reported in this thesis, other phase-domain frequency-dependent transmission line models have also recently been proposed and published as well [38 - 40]. 20 CHAPTER 3 Modelling of Single-Phase Nonuniform Transmission Lines for Representation of Transmission Towers 3.1 Overview A line model with exponentially varying parameters is introduced to represent nonuniform transmission lines in this Chapter. From the travelling wave standpoint, the trans- mission towers are viewed as nonuniform transmission line structures. When the line parameters are assumed to vary exponentially, a set of two-port equations can be formed in the frequency domain, which contain frequency-dependent functions. These functions are then synthesized with rational functions of the minimum-phase-shift type. Utilizing a fast recursive convolution technique, the time-domain equations of the proposed model reduce to a form similar to those in Bergeron's method. Thus, the model is compatible with general electromagnetic transients programs such as the EMTP [1]. In the early stage of this thesis work, a tapered line model [22] was developed for the representation of transmission towers. It had shortcomings, however, inasmuch as it was only a high frequency approximation. Subsequently, improvements were made with the exponential variation assumption of [23]. With this assumption, the non-uniformity of the line parameters can effectively be included in the time domain simulation, as described in the following sections. 21 z c ( x ) = Z e^ * 0 -> length = 1 <- Figure 3.1 - Single-phase exponential line 3.2 Exponential transmission line Notations: The uppercase represents frequency domain quantities, while the lowercase indicates their time domain correspondents. 3.2a Equations in the frequency domain The basic equations of a nonuniform transmission line, expressed in the frequency domain, are (3.1a) (3.1b) V and I are the voltage and current phasors, and Z(x) and Y(x) are the space-dependent per-unit-length series impedance and shunt admittance, respectively. Following the procedure of [23], Eq. (3.1a) is differentiated again, dV dx 2 2 y ( . d l , dZ(x) dx dx v Substituting (3.1a) and (3.1b) into (3.2) gives (3.2) 22 For the exponential line shown in Figure 3.1, with losses ignored, Z(x) and Y(x) are Z(x) = j(oL(x) = j®L e , (3.4a) qx 0 Y(x)=j(oC(x)=j(oC e-" , (3.4b) x 0 where L(x) and C(x) are the per-unit-length inductance and capacitance, respectively; L and 0 C are the values at x = 0. These parameters are related to the high-frequency approximation 0 of the line characteristic impedance, Z h i A x ) = where lm ° > =z - eqx (3 5) - 6 ) Z =fe. 0 Substituting Eq.'s (3.4a) and (3.4b) into Eq. (3.3) gives ^ where c = - ^ + ^ F = ' 0 ( 3 ^ is the wave speed. J LoC o Eq. (3.6) is a second-order differential equation with the roots x,.| Jfi) -( ) . , x-f-JS'-d) : 1 + J f (3.7) Then, V(x) = Cie^ x + C e^ , (3.8) x 2 and from Eq. (la), I(x) = -^[^C e^ x l + X2C2e nx The constants C, and C depend on the boundary conditions. 2 evaluating Eq.'s (3.8) and (3.9) at nodes k and m. At node m we have (3-9). They can be found by 23 xi V = C e^ + C e^', l m ] 1 I =-^=-[X C e^ + l m where Z i l (3.10) 2 2 XCe% (3.11) x 2 2 =jwL e . qx m 0 Thus, from Eq. (3.10), Cx=\V -C e^ \e-^ . l m (3.12) 1 2 Substituting Eq. (3.12) into Eq. (3.11) and collecting terms, A.) V +Z I X\ — X C m 2 m m (3.13) 2 Then C, is obtained from Eq. (3.12), X Vm 4~ Z I [_ X — X\ 2 m (3.14) m 2 Evaluating Eq.'s (3.8) and (3.9) at node k leads to v = c,+c , k (3.15) 2 1 h = -—[X\C\ + X C ], 2 (3.16) 2 where Zk = j®L . 0 Multiplying Eq. (3.16) with t — and adding it to Eq. (3.15), the voltages and currents at K 2 both ends can be related as [Vk + Z ((o)I ]A = V +Z ((o)I , ck where A = e Xxl k m is the propagation function. Z c k cm and Z (3.17) m c m are the characteristic impedances at both ends of the line, which can be expressed as j<nL e ' Z ((o)=X q Z k((a) = ^~ J c d an 0 cm (3.18) 2 By letting co go to infinity in Eq. (3.7) and inserting the result into Eq. (3.18), the characteristic impedances become the high-frequency approximation of Eq. (3.5). In [22], this high-frequency approximation was used. 24 3.2b Equations in the time domain Reversing the current direction at node m to make it flow into the line, I mk = -I , m Eq. (3.17) in the time domain becomes M O + z (t) * i (t)] * a(t) = v (f)-z {t) * i (i), ck km m cm mk (3.19) where the symbol * denotes convolutions. If the characteristic impedances Z c k and Z c m , together with the propagation function A , are synthesized with rational functions, their corresponding expressions in the time domain will become simple sums of exponential functions. The function a(t) will also have a time delay, which approximately equals the time it takes for the fastest frequency component to travel along the line. Accordingly, the convolutions in Eq. (3.19) can be evaluated with a fast recursive algorithm [34]. Thus, the voltage at node m can be expressed in a simpler form as v (t) = z i k(i) m where z meq meq m + v c (i) + Vh ro (f), hR m P m (3.20) is a constant. The last two terms on the right- hand side are evaluated from the known values of previous time steps. The first term comes from the RC-network which approximates the characteristic impedance, and the second term comes from the propagation of conditions at the remote end k. Combining these two terms into a single history voltage source e . (t), Eq. (3.20) becomes h m v (t) = z i (f) m meq mk + e -m(f). h (3.21) The same procedure is used to evaluate the voltage at node k, but with the wave direction from node m to node k instead. Factor q has opposite sign now, and the roots in Eq. (3.7) must therefore be re-evaluated. The characteristic impedances and propagation function of 25 k ^keq -> I W v ~ km e h-k' ^meq m V \ A A — — i <— + h -m m Figure 3.2 - Equivalent circuit for the exponential line the line are also recalculated, using these new roots. Carrying out the necessary steps, the voltage at node k can then be expressed in the same form as Eq. (3.21), v (t) = z i (t) + e - (t). k keg km h (3.22) k Eq.'s (3.21) and (3.22) are very similar to those in Bergeron's method [1]. They are compatible and can easily be interfaced with the EMTP, with the equivalent circuit of Figure 3.2. 3.3 Synthesis of exponential line functions Recursive convolutions leading to Eq. (3.21) and (3.22) are possible only if the line characteristic impedances and propagation functions can be synthesized successfully. thesis results are discussed in this section. The line parameters are synthesized using rational functions of real negative poles and zeros. high-frequency characteristic impedances Z h i g h k and Z A 50 m long exponential line with h i g h m of 220 and 150 Q at node k and node m, respectively, is used for demonstration purposes. 3.3a Travelling from node k to node m The factor q for the example is evaluated as q 1 Zhighm I Z h i g h k 1 , 50 The syn- 150 220 = -0.0076 26 Figures 3.3 and 3.4 show the magnitude and phase angle of Z (co)'. Note that it is only ck necessary to calculate and synthesize Z (co) because Z (co) can be found by multiplying ck cm Z (co) with e . Figures 3.5 and 3.6 show the magnitude and phase angle of the propagation ql ck function A . The fitting accuracy is quite reasonable in these figures, and the magnitude and phase of the synthesized functions trace those of the exact ones very closely. The fitting accuracy can be further checked by comparing the results directly in the time domain. Figure 3.7 shows the exact and synthesized propagation function a(t) obtained from taking the inverse Fourier transform of A(oo), using the Fourier transformation program of [55]. To filter out high frequencies, the Lanczos filter of the form ^ ( s n 7ia)/ ^ ° ax) ( m w a s a ppii j t ec Q the (710)/© ax) m synthesized and exact functions, where co is the truncated frequency and is equal 27tl0 for 8 max this case. The two curves in Figure 3.7 are practically identical, which illustrates that the fitted a(t) is an excellent representation of the exact one. The propagation function a(f) has a physical meaning [36, p. 4-95]: If the line were to be energized with a voltage from node k which has an amplitude of 1.0 at all frequencies, through an impedance matching the line characteristic surge impedance Z (co), and with ck node m open, the voltage V at the receiving end would be the propagation function A(co). m This, in the time domain, is interpreted as applying an impulse voltage (infinite magnitude with an area of 1.0) at node k through the characteristic impedance to the line. The voltage arriving at the open end of the line will then be the propagation function a(t). Note that the line under consideration is 50 m long, which means that it should take approximately 16.67 ps for the impulse to reach the other end of the line. This is indeed seen in Figure 3.7. 'it was found experimentally that the fitting of Z (co) is better in the low frequency region when a small conductance G is ck included in Eq. (3.4b). In that case, the term (co/c) in Eq.'s (3.6) and (3.7) becomes - Z Y . 2 27 x 1 0 • i • 1— i Exact function 5 • k \ - a 300 o - o 250 3 • 400 350 4 1 o 200 1 - 150 2 1 00 io 10 4 1 5 > 10 6 10 7 1 o 8 E n l a r g e m e n t n - 2 1 1 0 1 0 0 1 0 1 0 1 0 Frequency (Hz) Figure 3.3 - Magnitude o f Z c k in Q 5 0 Exact function Synthesized function 0 5 0 - 1 0 0 10 2 0 10 10 10 Frequency (Hz) Figure 3.4 - Phase o f Z c k in degrees 1 0 1 0 .28 0.85 Exact function Synthesized function 0 .8 0.75 0 .7 0.65 - 2 1 0 10 1 0 1 0 1 0 1 0 Frequency (Hz) Figure 3.5 - Magnitude o f A 2 0 0 - 2 0 Exact function E n I a rg em Synthesized function - 4 0 - 6 0 - 8 0 1 0 0 1 2 0 1 0 1 0 10 1 0 Frequency (Hz) Figure 3.6 - Phase o f A in radians 1 0 1 0 29 . 2l . —. 1 5 1 6 1 7 .—- .—. 1 8 1 9 1 2 0 Time (us) Figure 3.7 - Propagation function a(t) 3.3b Travelling from node m to node k As mentioned previously, when the wave travels in reverse direction, a different factor q must be used in the equations, which results in different characteristic impedances and propagation function. The prime symbol will be used to identify these functions, namely Z ' , Z ' , and A'. They are used in Eq. (3.22) to evaluate the voltage at node k. The factor ck c m q in this case is the negative value of what it was for the wave travelling from k to m, namely l l n |^!L = 0.00766. Figures 3.8 and 3.9 show the magnitude and phase angle of Z ' . The function Z ' c m lated to Z ' is re- with a multiplication factor of e . On the other hand, A' can be evaluated diql c m c k rectly from the propagation function found for the other direction. 2 0 0 Exact function Synthesized function 1 5 0. 1 0 0 5 0 0 1 0 10 10 10 1 0 1 0 Frequency (Hz) Figure 3 . 8 - Magnitude o f Z ' c m in Q 1 0 0 8 0 Exact function Synthesized function 6 0 4 0 2 0 0 1 - 2 0 1 0 10 10 Frequency (Hz) Figure 3 . 9 - Phase o f Z ' c m in degrees 1 0 1. 0 31 Z h (3.23) ighm It can be seen that A' differs from A only by a constant factor. Thus, once A is synthesized, the approximate function for A' is also known! The good agreement between the fitted and the exact A functions applies to A' as well. 3.4 Reflection factors To gain confidence in the proposed line model, it is useful to calculate its input reflection factor and compare it against previously published results. Eq.'s (3.8) and (3.9) are the starting point for this comparison. The impedance looking into the line at any point x, Z (x), can be computed by simply taking the voltage-to-current ratio from these two equain tions. Therefore, the impedance looking into node k, at the sending end, is evaluated with x = 0 as Z,„(0) = Ci+C -z 2 k _ X\C\ + XC 2 2 Using the constants from Eq.'s (3.13) and (3.14), (A.2 V + Z I )e-^ 1 Z (0) = Z in m k m m kifri V + Z I )e-* 1 m m m Z (x)= c -(XxV + m -Xi(X V + 2 m Z I )e-^ 1 m m Z I )e m m z eqx 0 length = 1 Figure 3.10 - Single-phase exponential line terminated at node m with Z L 32 Letting the exponential line be terminated at node m with a load impedance Z as shown L in Figure 3.10, then V = Z I , and.Z (0) becomes m L m in D e' ' x 1 Z (0) = Z in x k _X D e- ' k2 2 2 -D e~ ' X2 2 (3.24) -XtDie-*- '. 1 where Di=k Z 2 +Z L m and D =X Z 2 X L +Z . (3.25) m Note that Z and Z are the per-unit-length series impedances evaluated at node k and node k m m, respectively. Figures 3.11 shows the magnitude of the input impedances evaluated from Eq. (3.24) for two different values of Z . It is interesting to note that this input impedance is exactly the L same as the load impedance Z , at frequencies less than about 10 Hz. This is because the 5 L propagation function is purely real below this cut-off frequency, which presents only attenuation in the electromagnetic field without any progression (no travelling wave). Under this condition, the load impedance Z is seen as if it were connected right at the sending end of L the line. When travelling waves enter node k, the reflection coefficient seen at this node is defined as . ™=tf^i? <> 326 with the assumption that a uniform line of characteristic impedance Z is connected to the 0 sending end at node k. In [29], the reflection coefficient is approximated as 33 where (3 = Eq. (3.27) was derived from afirst-orderapproximation for the reflection coefficient between two infinitesimally short sections. It is therefore approximate and not as accurate as Eq. (3.26). In addition, Eq. (3.27) was established only for the case when the load impedance Z and the high-frequency characteristic impedance at the receiving end m are nearly L equal to each other [29]. impedance Z = Z L h i g h m For the exponential line considered in section 3.3, with the load = 150 Q, Figure 3.12 shows the magnitudes of the reflection coeffi- cients evaluated from Eq.'s (3.26) and (3.27). Both methods produce close answers for the line under consideration. The differences become larger when the Z h j g h k /Z h i g h m ratio increases to greater than 2. If Z is set to 100 Q instead, Eq. (3.27) is no longer valid, as shown in FigL ure 3.13. 3 5 0 8 0 10 1 0 . 10 Frequency (Hz) Figure 3.11 - Magnitude o f Z i n in Q from E q . (3.24) 10 34 0 . 2 10 10 10 10 10 Frequency (Hz) Figure 3.13 - Magnitude of Reflection factor for Z = 100 Q L 35 3.5 E M T P simulations To show the behaviour of the proposed nonuniform line model, EMTP simulations of the line responses under open- and short- circuit conditions were performed. The results are shown in Figures 3.14 and 3.15. For testing purposes, the equations of the exponential line were interfaced with the U B C EMTP version MicroTran through subroutine C O N N E C . The interface is described in Appendix C. Additional time-domain simulations for comparisons with field measurements will be shown in Chapter 6. An exponential line can be approximated as a cascaded line model [29], where the line is divided into many sections of uniform transmission lines. Figures 3.14 and 3.15 show simulation results for open and short circuit tests for this approximation as well. The exponential line, with parameters described in section 3.3, is energized with a 1 p.u. step voltage at node k. Three models are used for comparison purposes: 1) the new model, 2) a cascaded model consisting of 50 short uniform lines, and 3) an ideal uniform line having a uniform characteristic impedance equal to the average value (185 Q). Note that the proposed new model agrees well with the cascaded model for the open circuit voltages. It clearly shows the tapered effect created by the continuous reflections along the line, in contrast to the square waves obtained from the uniform transmission line model. The voltage obtained from the new model has larger peaks, compared to those produced by the cascaded model, because the latter model does not have enough short sections to represent these peaks. The short-circuit currents agree more or less among the three models during the first 1.5 ps. The current from the new model then starts to deviate from the others because of the continuous reflections along the line. 36 3 . 5 Cascaded model 3 U n i f o r m model 2 . 5 2 1 . 5 1 0 . 5 0 -0.5 - 1 N e w model -1.5, 1 0~5 1.5 2 T 2 Time (u.s) Figure 3.14 - Receiving end open-circuit voltage (p.u.) 0.07 0.0 Cascaded and U n i f o r m models 6. 0.05. 0.04. N e w model 0.03. 0.02. 0.01. 0 o: 5 —\— 1 r^T" Time (us) Figure 3.15 - Receiving end short-circuit current (p.u.) 37 CHAPTER 4 Frequency-Dependent Transformation Matrices for Untransposed Transmission Lines using Newton-Raphson Method 4.1 Overview The development of a complete frequency-dependent overhead transmission line model in the time domain is one of the major objectives of this thesis. In the frequency-dependent overhead transmission line model, presently in the EMTP, constant and real transformation matrices are used for conversion between modal and phase domain quantities. As men- tioned in section 2.4, constant transformation matrices are only accurate enough for single flat circuit lines. For the case of double- or more circuit, lines, the accuracy decreases noticeably, especially in studies which cover wide frequency ranges [33]. Lightning is a phe- nomenon which, in most cases, involves a wide range of frequencies. For example, in the simulation of the shielding-failure flashover and back flashover events, the 60 Hz steadystate as well as the transient voltages are needed. Thus, a wideband frequency-dependent line model is needed for such situations. One such model is the full frequency-dependent modal-domain model of [35]. It includes the complete frequency-dependent nature of the transmission line and synthesizes not only the modal characteristics but the elements of the transformation matrices as well, with rational functions. To achieve this, the elements of these matrices must be continuously and smoothly evaluated over a wide frequency range 38 (1 Hz - 1 MHz). A stable eigenvector/eigenvalue routine based on the Jacobi algorithm was developed and used in [35]. The smoothness and minimum-phase requirements for the transformation matrix elements were satisfied, but only for the case of underground cables. For multi-circuit overhead lines, however, these elements could not be fitted with rational functions of minimum-phase shift. In addition, they did not always vary as smoothly as one would wish [45]. Realizing this limitation, a diagonalization routine was developed in this thesis project, together with an effective scaling method, which produces well-behaved transformation matrices. The time-domain implementation procedure can then follow the same formulation as described in [35]. The diagonalization method is based on the New- ton-Raphson iteration technique. It was first proposed by L . M . Wedepohl [48], and initially programmed by G.D. Irwin at the Manitoba H V D C Research Centre. The method was then additionally improved and implemented successfully with the U B C Line Constants Program by the author [37]. This algorithm produces well-behaved transformation matrices for overhead transmission lines. Although some elements of the transformation matrices may present non-minimum phase behaviour, they could still be synthesized with rational functions which have real negative poles and zeros, with reasonable engineering accuracy. 4.2 Equations for multiphase transmission lines Notations: Boldface represents matrix quantities, whereas uppercase depicts frequencydomain quantities. Considering a multiphase transmission line system as shown in Figure 4.1, the basic transmission line equations in the frequency domain are: 39 Node k Ik Node m Im legnth = l > n • > > Phase #1 V Phase #2 •n o —> Phase # n Figure 4.1 - Multiphase transmission line system - f = ZI, - ^ = YV, ax (4.1.). (4.1b) where Z and Y are the per-unit-length series impedance and shunt admittance matrices, respectively. V and I are voltage and current vectors (or column matrices). dV Differentiating Eq. (4.1b) again and replacing — with Eq. (4.1a), we get for the line ax currents, V = [YZ]I. ax (4.2) It is well known that Eq. (4.2) is solved by transforming the n coupled equations into n decoupled equations. To achieve this, a modal transformation matrix T „ which diagonalizes the Y Z product, is found from the solution of the eigenproblem, T " Y Z T i = A,, 1 I (4.3) where A, is the diagonal eigenvalue matrix. The matrix T, varies with frequency since Y and Z are frequency dependent. 40 In order to understand fully the frequency-dependent characteristics of the matrix T „ it is best to examine the frequency-dependent nature of Y and Z first. The matrix Y has the following form [36, pp. 4-16 to 4-17], Y = G+yooP-', (4.4) where co = 2nf is the angular frequency, G is the per-unit-length conductance matrix, and P is the constant Maxwell potential-coefficient matrix, having each element defined as s = permittivity of free-space, 0 Djj = distance from conductor i to image of conductor j , dy = direct distance form conductor i to conductor j , or conductor radius when j = i. The conductance matrix G is typically represented as a diagonal matrix. G is small compared to jcoP" , and is therefore often ignored. In [34] and [36, p. 4-99], values of 0.03 1 uS/km were used. The matrix Z has the form of Z=y'cop s P + Ze + Z , 0 where 0 (4.6) c is the permeability of free-space. Subscripts c and e refer to the internal impedance of the conductor and earth impedance corrections, respectively. Note that Z is a diagonal c matrix. For solid conductors, its elements can be evaluated with a simple but sufficiently accurate formula in the frequency range of interest here [49], -f PO Zcii = [ ^ - j I (i —— r icothi — 0.777AI +0.356 ^ , c ' \2d \J V d\ J n e e A (4.7) 41 where: p, = conductor resistivity, r = conductor radius, del = I , conductor complex depth of penetration, p = conductor permeability. At very low frequencies, d . approaches infinity, and Z u = -—r, which is the correct d.c. c resistance. At very high frequencies, d el Pi approaches zero, and Z u — > -——, 2nrd ] c which is the e correct skin effect impedance at high frequencies. The formula for evaluating the earth correction impedance matrix Z was developed by e Carson [50]. This formula has been used by power engineers for many years. Deri etal evaluated the elements of this matrix with a simpler formula [51], which was originally proposed by Dubanton [52]. It has the following form [36, pp. 4-4 to 4-5] and [41], 2 * - ^ * $ . lit JJii (4.8) where D'j = ^(y,+y, + 2af ) + (*,• -Xj)' 2 e2 The term d i = / ~— e is the earth complex depth of penetration, x and y are the hori- zontal and vertical co-ordinates of the conductors, respectively, and p is the earth resistivity. 2 Eq. (4.8) is accurate enough for engineering purposes, and the results agree very closely with those obtained from Carson's integral formula. The maximum error is only a few percent in the kHz range. This is of very little consequence in transient calculations where the entire frequency spectrum contributes to the results anyhow. Thus, Eq. (4.8) is preferred in this thesis because of its simplicity. It should be mentioned that the transformation matrices 42 evaluated in this thesis, using Carson's formula or Eq. (4.8) in the Y Z product, were practically identical. 4.3 Eigenproblem solution with Newton-Raphson (NR) method For convenience, the following additional variables are introduced: i) Sis the Y Z product, ii) T I ( k ) is the column k of matrix T , , iii) U is the unit matrix. Then, equation (4.3) can be expressed as ST, =m . (4.9) Eq. (4.9) contains n equations for an n-phase system, because each of the n eigenvalues 2 X is related to n equations. Considering one set of n equations, which belongs to a particuk lar eigenvalue X [41], i.e.: k (S-X U)T k m = 0. . (4.10) Eq. (4.10) is a set of homogenous simultaneous equations. It may be shown that the solution is trivial unless the determinant of (S - A,*U) is zero. Thus, the matrix S can be diagonalized by equating this determinant to zero. This leads to a polynomial in X of order n, yielding n values of X. Sometimes one or more values of X may be equal. Once the values of X are known, the elements of T I ( k ) may be found from Eq. (4.10). There are various software packages available for determining the eigenvalues and ei-. genvectors. However, for the modal-domain transmission line model of [35], the eigenva- lues and eigenvectors must be found at many frequencies. The elements of the eigenvectors 43 have to vary smoothly as a function of frequency throughout the entire frequency range. This is very important because rational approximations for these elements are required. If a normal algorithm is used, the eigenvalues will switch places at certain frequencies, which creates jumps in the eigenvectors. A complicated "tracking procedure" was used in [35] to avoid jumps. The N R algorithm proposed here overcomes this limitation. It solves Eq. (4.10) iteratively for each value of A, together with each set of T , k tions has "n +1" unknowns. (k) . This system of non-linear equa- The number of variables could be reduced by one, by specifying an arbitrary value for one of the elements of T I(k) . For example, the element with the largest value at the initial frequency could be chosen and set equal to 1.0 for all frequencies. However, this could result in overflow or excessive errors, because the value of this element may sometimes become zero or very small at some frequencies. A safer way is to introduce a constraint equation. One such constraint is to specify that the sum of the squares of the elements of the eigenvector be unity. This constraint method has been chosen in this thesis. With this constraint, the Jacobian matrix is well behaved, and the rational approximations of the resulting smooth transformation matrices are found with ease. For example, for a three-phase line, one set of 3 equations belonging to A,, and column 1 of the matrix T , from Eq. (4.10) are (5,i-Xi)T =0, (4.11) 5217-/1,+(522-^)7/21+523^1=0, (4.12) 53ir/„+532r/2,+(533 -Al)T/31 =0. (4.13) m +SuT m +S T n Dl 44 There are four unknowns but only three equations. The extra equation is the constraint equation, which specifies that the sum of the squares of the elements of T I ( 1 ) be 1.0, T +T +T -l=0. 2 2 m (4.14) 2 m m The Jacobian matrix for the above four equations then has the form, (S]\-X\) J = Sn S\3 -T/ii 521 (^22-A-l) 523 -721 531 532 27/n 27/2i (533 - ^ l) (4.15) -T31 27/3i 0 Let x represent the column vector of the four unknowns Tn\,Ta\,Tn\ and V Then x in iteration step n is found from x (n) _ ( " - l ) . x j( («-i)) ' F ( x " - ) ( 1) X (4.16) where J(x)" is the inverse of the Jacobian matrix and F(x) represents the left-hand sides of 1 equations (4.11)- (4.14), which eventually become zero as the iterations converge. The iteration procedure starts from an initial guess x . (0) As starting values of x, the solution values from the previous frequency are used. The NR method generally converges very fast, provided that reasonable starting values of x are known and that J(x)" exists. After convergence is reached, A, and column 1 of the matrix T, 1 are found. The same procedure is then followed for X. and column 2 of T , , and then for X 2 3 and column 3 of T , : It should be noted that this N R method is not self-starting. Another method must be used to supply the solution for the initial frequency first. In the program written for this thesis, the diagonalization routine based on the power method provided by L . M . Wedepohl [41] is used to start at the initial frequency. In addition, in order to avoid floating-point overflow 45 and to ensure continuity of the modal parameters at the high frequency region, the Y Z product is scaled before the N R diagonalization algorithm is called. The scaling procedure consists of dividing each element of this product by a constant -(co p s ), and subtracting a unit 2 0 0 matrix (1.0 for diagonal and 0.0 for off-diagonal elements) from the resulting matrix. As seen in Appendix A , this scaling procedure does not alter the eigenvectors of the original Y Z product. This scaling procedure was suggested in [54], and is further explained in [36, pp. 4-52 to 4-53]. 4.4 Modal parameters for multiphase transmission lines from Newton-Raphson method Modal parameters were found with the N R method for different line configurations, ranging from single- and double- circuit lines for horizontal and vertical constructions to uncommon circuit lines with the phase conductors arranged in various asymmetric settings, with or without ground wires. A n example of an uncommon double-circuit line is the case where the left three phases are vertical while the right ones are horizontal. Another case with strong asymmetry, which is common in practice is the case where two or more circuit lines of different voltage ratings share the same right-of-way and may be located one underneath the other. The N R method had no difficulty in diagonalizing such uncommon line configurations, and the transformation matrices were all smooth and well behaved. The fitting routine used in these tests, referred to as FIT-S, is a modified version of the one described in [35]. This fitting routine was developed by J.R Marti [34]. The main modification involved the calculation of the time delay constant for the rational approximation of the modal propagation function. Since the angle of this function can be found 46 directly, and since it is smooth and continuous, only one time delay at the fastest frequency was used, whereas the average time delay over a frequency range was used in [35, 45]. The synthesis and fitting criteria for all modal functions are explained clearly in references [34] and [35]. They will not be repeated here. It is significant to mention that the N R method produces modal parameters which are smooth and vary only slowly as a function of frequency, because of the inherent eigenvalue-eigenvector tracking nature of the N R algorithm. This is the reason why syntheses using minimum-phase-shift rational functions have been successful in this thesis (FIT-S has not failed in fitting the parameters produced by the N R method for all the line configurations tested in this thesis). 4.4a M o d a l transformation T , matrix There are n modes in an n-conductor transmission line. In the simulation of lightning tower/ground wire surges and back flashover events, ground wires are included as conductors, whereas in shielding failure and shielding-failure flashover events they may be eliminated from the Z and Y matrices by Gaussian elimination. For the double-circuit line of Figure 4.2, assuming that the ground wires have uniformly zero potential and can thus be eliminated, there are six modes describing this line. Figures 4.3 shows the magnitude and phase angle of the elements of column 1 of the 6 x 6 transformation matrix. Note that there are actually twelve curves superimposed in each of these, plots. It is seen that the fitting functions match the exact functions closely. The elements in other columns also match quite reasonably, but some differences in the angles showed up in the low frequency region. 47 , 11.1 m A . G r o u n d wire: 1.75 c m d i a m e t e r 7fT _\/_ 20.8 m . 1 19.0 m r O c . length = 100 km n 19.0 m 6 4 Phase conductor: 4 X 3.84 c m d i a m e t e r at 50 cm spacing 12.7 m 2 O 20.8 m 12.7 m O 5 o- o 3 O o 6 ! I • i i o—o 50 cm 26.0 m 100 Ohm-m ground resistivity Figure 4.2 - Vertical double-circuit line These differences, however, have very little bearing on the overall transient solution, as proven in the open/short -circuit tests in section 4.5. Because of symmetry, the elements of each column of the matrix T , occur either in equal or in equal and opposite pairs. The phase conductors are arranged such that elements 1 and 4, 2 and 5, and 3 and 6 of each column are equal in magnitude. It is shown in Appendix A that the elements of the matrix T , become real and approach asymptotic values at the high and low frequency ends. This characteristic is seen in Figure 4.3. 4.4b Modal characteristic admittance matrix Y,C M The modal shunt admittance and series impedance matrices are defined as Y M = T7'YT7 r and Z M = T[ZTi , (4.17) where T," is the inverse of the transposed matrix T, . T T It follows from Eq. (4.9) that Y M Z M = A,. (4.18) 48 It has been proved that both Y mittance Y C M and Z M M are diagonal [53]. The modal characteristic ad- is also a diagonal matrix and can be found from the scalar relationship, (4.19) YcM-k where Y C M . is the diagonal element k of the Y k C M matrix. Figure 4.4 shows the six elements of the modal characteristic admittance matrix and their corresponding synthesized functions for the line in Figure 4.2. Note again that the two fitted and exact curves are superimposed in the figure. 4.4c Modal propagation matrix A M The modal propagation matrix A M is also diagonal. Its elements have the form of AM-k = e-J * = e- e-W, 1 1 (4.20) al where X is the k-th eigenvalue and 1 is the,length of the line. The term e" is the attenuation al k factor and e" the phase shift. Figure 4.5 shows the superimposed (exact and synthesized) jpl curves of these matrix elements for the line in Figure 4.2. There are only five distinct el- ements shown in each of the plots since two of the eigenvalues are very close or identical for this line. The phase angles, without time-delay phase shifts, are shown as well. Again, the fitting is very good. In order to ensure continuity in the travel time delays, it is important to obtain the phase angle of each element A . directly as M k arg{A -k} = Im{JTk}l. M (4.21) 49 With this equation, the angle is smooth and continuous, as shown in Figure 4.6. were evaluated as the phase angle of the function 1 If it , using the intrinsic function atan in F O R T R A N , it would show highly oscillatory behaviour in the high frequency region, as shown in Figure 4.7. It would be difficult to evaluate the time delays needed for the fitted functions with such a behaviour in the phase angles. This aspect will be elaborated further in Chapter 5. FIT-S usually allocated fewer poles for the elements of Y C M and T , than for those of A . On the average, twelve poles were used for the former and 28 for the latter. M 50 Exact function 0 . 8 Synthesized function 3 &6 0 . 6 <D 3 CD ro 1 &4 Maximum error < 1.0 % 0 1 0 10 1 0 1 .0 Frequency ( H z ) 2 0 1 0 in T3 D) 0 E c ro (U (/> ro 1 0 2 0 1 0 1 0 1 0 Frequency (Hz) Figure 4.3 - Elements o f column 1 o f the matrix T , 1 0 51 , X b , 10 , 1— —. # 0 10 2 10 Frequency (Hz) Figure 4.4 - Elements o f the matrix Y , 10 52 10 10 10 Frequency (Hz) Figure 4.5 - Elements of the matrix A , 10 53 x 1 0 m c ro ~o ro _0) TO C ro w ro 1 0 1 0 10 10 1 0 1 0 Frequency (Hz) Figure 4.6 - Phase angle of element A _ , calculated from E q . (4.21) M ro a> a) c ro <u w ro o k 1 0 1 0 10 10 1 0 Frequency (Hz) Figure 4.7 - Phase angle o f element A . , calculated with atan function M 1 0 54 4.5 Open/Short -circuit tests in the frequency domain An efficient and simple way of comparing the solutions produced by the exact and approximate modal parameters is to calculate the line responses under unfavourable conditions, such as short- and open- circuit tests, in the frequency domain. In order to perform these tests, a column voltage matrix E is applied at the sending ends (node k) with the res ceiving ends (node m) either open or shorted. The single-phase open-circuit voltage and short-circuit current for each mode at the receiving ends are evaluated as follows [45], ^(«» =T ^ # ji , x / _/(co) = m 1 . (4-22) 2YcM-iA -iE' _i M —— s , (4.23) 1 ~Am-I where Y . j and A . j are the i-th elements of the modal characteristic admittance and propaCM CM gation function matrices, respectively. The prime symbol indicates that the voltages and currents are in the modal domain. The conversion between modal and phase domain for the voltage and current column matrices is made with the following relationships, I (co) = T,(co)I' (co) and V (co) = T7 (coVVUco). r m w m The exact line responses are calculated using the exact T , , A and Y M ( 4 C M . 4) 2 matrices, where- as the approximate responses are obtained using their respective approximations by rational functions shown in section 4.4. In the following, two forms of excitation will be used to test the double-circuit line of Figure 4.2. 55 4.5a Zero sequence excitation For this form of excitation the source voltages all have 1 p.u. magnitude and zero angle. Figures 4.8 and 4.9 show the magnitude of open- and short- circuit responses for all phases at the receiving end. Note that the approximate and exact responses agree very closely. The angles are not shown but they also agree closely. 4.5b Unsymmetric excitation The accuracy of the line responses with approximate modal parameters depends on the degree of imbalance in the voltages and currents. A n unfavourable case of voltage source imbalances is used to test the line considered here. This case was suggested in [33], with these values: E \ =E 4 = IZ0°, s S E =E s =0, and E& =E 6 = 1Z120 . 0 s2 s S Figures 4.10 and 4.11 show the magnitude of the responses for all phases at the receiving end. Note again that the approximate and the exact solutions match one another almost perfectly. These tests were repeated for many other line configurations. The degree of accuracy is generally very good for the cases of single- and double-, flat and vertical, circuit lines. For lines with strong asymmetry, the approximate and exact solutions still agree quite well when the excitation sources are unsymmetric or symmetric. For the zero sequence excitation, small differences appear at some frequencies in the medium frequency range. For all practical purposes, these small differences are still considered to be acceptable in the simulation of power system transients. Exact function Synthesized function - A. - 2 1 0 1 0 0 1 0 2 1 0 4 1 0 6 1 0 Frequency (Hz) 1 0 8 6 4 2 -r-~ 1 0 n r 1 0 1 0 1 0 1 0 1 Frequency (Hz) 1 0 1 0 Frequency (Hz) Figure 4.8 - Receiving end open-circuit voltages (p.u. magnitude) ( 57 Figure 4.9 - Receiving end short-circuit currents (p.u. magnitude) 7 0 6 0 Exact function Synthesized function 5 0 °3 4 co a) CO ro CL 0 3 0 2 0 1 0 1 0 1 o 1 o Ll 1 0 1 0 1 0 Frequency (Hz) 5 0 4 0 - LO CN w w ro 3 0 - 1 0 10 1 0 10 Frequency (Hz) 5 0 4 0 - CO oa CO CO (U CO ro 10 10 1 0 Frequency (Hz) Figure 4.10 - Receiving end open-circuit voltages (p.u. magnitude) 1 0 1 0 . 1 1 0 0 1 0 1 0 Frequency (Hz) 0 . 2 : : \ Figure 4.11 - Receiving end short-circuit currents (p.u. magnitude) 60 4.6 Validation of the Newton-Raphson algorithm with the power method An alternative method was used to validate the results obtained from the N R algorithm. This alternative method uses a diagonalization routine based on the power method. Since the eigenvalues will most likely switch places at certain frequencies, a correlation technique is used to keep track of the eigenvectors throughout the frequency range. This correlation technique is based on the fact that the eigenvectors belonging to the same set of eigenvalues are nearly orthogonal from one frequency to the next. This means that if an eigenvector at the present frequency is to be correlated with each and every eigenvector obtained at the previous frequency (the most recent frequency), it will be found that the strongest correlation factor will be produced by the two eigenvectors belonging to the same set of eigenvalues. The correlation procedure for tracking the eigenvectors can be summarized as followed: i) Obtain the eigenvector matrix T for the present frequency. I 2 ii) Transpose T I 2 and obtain its complex conjugate T iii) Find the matrix product T I 2 I 2 . T n , where the eigenvector matrix T„ is the one from the previous frequency. iv) In each row of the product in iii, identify the largest element. The row number of this element indicates the column number of the eigenvector from the previous frequency, where as its column number represents the column number of the eigenvector at the present frequency. For example, in the case of a 3 by 3 matrix, if eigenvalue switching has not occurred, the correlation product of iii would have the largest elements in the diagonal locations as indicated by the X's in Figure 4.12a, whereas if the eigenvalues of elements 1 and 3 have been switched, there will also be a switch between the eigenvectors of columns 1 and 3, 61 and the correlation matrix in iii would have the largest elements in the locations shown in Figure 4.12b. At every frequency, after the eigenvectors are obtained from the diagonalization routine, they are correlated with those obtained at the previous frequency. When switchings among eigenvectors are detected, they have to be rearranged in the order of the eigenvectors from the previous frequency. This guarantees their continuity throughout the frequency range. Tests confirmed that the modal transformation (eigenvector) matrices obtained from the NR method and from the above alternative method agree perfectly. Thus, the validity of the proposed N R method was proven. X X X X X (a) X (b) Figure 4.12- Location of largest elements of the correlation matrix 62 CHAPTER 5 Frequency-Dependent Overhead Transmission Line Model in the Phase Domain 5.1 Overview As an alternative to the modal-domain transmission line model described in reference [35] and in Chapter 4, the model presented here works directly in the phase domain. This approach has the advantage of not having to synthesize the transformation matrices, and also uses fewer numerical convolutions at each time step. Thus, it is computationally more efficient than the modal-domain models of untransposed lines. Direct phase-domain modelling has been proposed in [38] - [40]. The approach developed in this thesis differs from these references in that it synthesizes the frequency-dependent line functions directly in phase co-ordinates, using stable rational function approximations with real negative poles and zeros. This new wideband transmission line model has been tested extensively for various overhead lines with different and asymmetrical geometries. It is stable and provides answers which are accurate enough in both steadystate and transient conditions. 5.2 Phase-domain transmission line modelling Notations: Bold and normal faces represent matrix and scalar variables, whereas uppercase and lowercase are used for frequency-domain and time-domain quantities, respectively: 63 Considering again the multiphase transmission line system of Figure 4.1, the classical solution for Eq. 4.2 in the frequency domain is I = e- I +e l , rx x where T = J\Z (5.1) rx f b is the propagation function, and I and I stand for the vectors of forward f b and backward travelling currents, respectively. The voltage vectors of the line can be found as V* = - Y " ^ = Y c V * I / - <? h), 1 r (5.2) rx where the characteristic admittance Y is c Y = J(YZT Y l c (5.3) Pre-multiplying Eq. (5.2) with Y and adding the result to Eq. (5.1) gives c Y V , + I = 2e- I/. (5.4) rjc c x Applying boundary conditions at both ends, At node k, x = 0: Y V c k At node m, x = l: Y V m c + I = 2I/, (5.5) k -I = 2e- I . (5.6) rl m f In Eq. (5.6), I is the current flowing into the line, or !(*=[) = - I . m m Substituting Eq. (5.5) into Eq. (5.6), Y V -I =A(Y V c m m c + I ), k k (5.7) where A = e~ . (5.8) rl Eq. (5.7) can be rewritten as I =YcV -A(Y V m m c k + I ). (5.9) + I ). (5.10) u Analogously for the currents at node k, Ik = Y V - A ( Y c k c V m m 64 Y and A in Eq.'s (5.3) and (5.8) are the line phase-domain characteristic admittance and c wave propagation matrices, respectively. They can be evaluated using modal parameters obtained from the eigenproblem of Eq. (4.3) [41, 35], Yc = T , Y c M T , , (5.11) A = T,A T7'. (5.12) r M Note that the modal characteristic admittance matrix Y trix A M C M and the modal propagation ma- are calculated from Eq.'s (4.19) and (4.20). Transforming into the time domain, Eq.'s (5.9) and (5.10) become i«(0 = y (0*v«(0-a(0*f*(0. (5-13) i*(0 = y (0 * v*(f) - a(0 * f (t), (5.14) f*(0 = y (0*v*(0 + u(0, (5-15) f«(0 = y (0 * v „ ( 0 + i (t), (5.16) c C m where c c m and the symbol '*' denotes matrix-vector convolutions. If the elements of Y and A are synthesized using rational function approximations, then c the corresponding time-domain functions y (t) and a(t) become matrices whose elements c contain simple sums of exponential functions [34, 35]. The convolutions in Eq.'s (5.13) and (5.14) can then be evaluated using fast recursive methods [30, 34, 35]. In terms of the synthesis functions, the nodal equations at both ends of the transmission line can be expressed in the more convenient form (Appendix B), im(f)=yi V (t)-ihta-m(t), q m i*(f) = YeqVk(t) ~ his,-k(t), . (5.17) (5.18) 65 k m Figure 5.1 - Multiphase transmission line equivalent circuit in the time domain where y eq is a real, constant symmetric matrix. The second term on the right-hand side of Eq.'s (5.17) and (5.18) is the history current source evaluated from the known values of the variables at previous time steps. These equations are represented with the equivalent circuit of Figure 5.1, which can easily be incorporated into the EMTP. 5.3 Synthesis of the propagation and characteristic admittance matrices The N R diagonalization algorithm described in Chapter 4 was used to evaluate the phase-domain characteristic admittance and propagation matrices as functions of frequency. The syntheses were then performed using the routine FIT-S. 5.3a Synthesis of the characteristic admittance matrix Y c The characteristic admittance matrix Y in the phase domain is a square symmetric mac trix. Therefore, only n(n+l)/2 elements need to be synthesized. Experience in synthesizing this matrix shows that all of its elements can be fitted using minimum-phase-shift functions, that is, using only real negative poles and zeroes. The only detail is that the off-diagonal elements need to be multiplied by -1.0 before the synthesis is performed. The reason is that, 66 with this multiplication, the value of the modified function asymptotically becomes a real constant with zero phase angle at the high frequency end (Figure 5.2), which allows FIT-S to allocate the constant R (the function's high frequency limit), the poles and the zeros with 0 better accuracy. Figure 5.2 shows the magnitude and. phase angle of the elements in column 1 of Y for the vertical double-circuit line shown in Figure 4.2. Because negative c functions of the off-diagonal elements were synthesized, the phase angle of the off-diagonal elements shown there differs by 180 degrees from the original functions. Note that the synthesized and exact functions match very closely, with errors typically under 2% within the required frequency range (1 Hz - 1 MHz). If the off-diagonal elements were fitted directly without the multiplication with the factor -1.0, noticeable differences would show up between the exact and synthesized functions at the low frequency region. 5.3b Synthesis of the propagation matrix A The intention is to synthesize all the elements of the propagation matrix A with minimum-phase rational functions, and associate a different time delay with each of them. To achieve this, some function manipulations are required. First, the modal diagonal matrix A M is evaluated from Eq. (4.20). Secondly, a phase shift is then extracted to give PM = ^ A M ; (5.19) where x is the travel time of the line based on the speed of light. Finally, a triple matrix product is calculated, P = TiP T7'. M (5.20) 67 The resulting matrix P is then synthesized using rational functions. Each of the syn- thesized elements contains an additional (small) phase shift which corresponds to an additional time delay in the time domain. As in the case of the matrix Y , the off-diagonal c elements of the matrix P must be multiplied by -1.0 before the synthesis is performed. Figure 5.3 shows the fitting of the magnitude and phase angle for the elements of column 1 of the matrix P, for the double-circuit line of Figure 4.2. It is seen that the exact and fitted magnitudes agree closely, and the phase angles are in good agreement as well. Because.of the phase shift extraction in Eq. (5.19), the angles of the elements of the matrix P are smooth throughout the required frequency range (Figure 5.3), and the additional time delay can be accurately determined. The approximation for the propagation matrix A is found as &fitted = e- ?fittedJW (5.21) Note that each element of the synthesized matrix A has a single time delay. This time delay is the sum of the delay x of the speed of light and the delay from the fitted matrix P. The multiplication by -1.0 of the off-diagonal terms of the matrix P affects very much the synthesized functions. Without performing this, the phase angles of these elements would be 180 degrees from the minimum-phase behaviour, as shown in Figure 5.4a. With this, it is extremely difficult to find the correct time delays for the fitted functions. There may be a physical meaning which justifies the multiplication of the off-diagonal elements with the factor -1.0: by observing the propagation functions in the time domain, as shown in Figure 5.5, one can see that when one phase is energized with an impulse, the other phases have negative surges through coupling. As explained in section 5.3c, the impulse 68 waveforms arriving at the other end of the line correspond to the elements of a column of the propagation matrix a(t). If the matrix A were to be calculated directly from the triple product in Eq. (5.12), without first extracting the phase shift e" , its elements would show highly oscillatory angles j(0T (between +180 and -180 degrees) in the high frequency region, as illustrated in Figure 5.4b. This would make it nearly impossible to relate the fitted function to a single time delay. For example, if A,,(co) were to be synthesized with FIT-S, a minimum-phase function, say F (co), would be obtained. The two functions are related to each other as H ^ii(ro) = e-- ' Fn(co), / fflw 9—0 where imax = co respectively. x max ^' a n c * w n e r e a n d @F a r e the phase angles of A,,(co) and F,,(co), is the.time delay calculated at the fastest frequency co , within the max considered frequency range. In order to evaluate this time delay with acceptable engineering accuracy, the phase angle of the A,, element has to be smooth and continuous throughout the entire frequency range. This condition would not be met if the phase angles were obtained directly from Eq. (5.12), as shown in Figure 5.4b. 69 Figure 5.2 - Elements of column 1 o f matrix Y c for a double-circuit line 70 10 10 10 10 10 10 10 Frequency (Hz) 10 10 10 Frequency (Hz) Figure 5.3 - Elements of column 1 o f matrix P 71 2 0 0 1 0 0 0 1 0 0 2 0 0 L— 1 0 z 10 10 1 0 1 0 Frequency (Hz) Figure 5.4a - Phase angle (degrees) o f elements P , and P , before multiplying with -1.0 2 3 2 0 0 10 0k 0 - 1 0 01 2 0 0 1 0 10 10 1 0 Frequency (Hz) Figure 5.4b - Phase angle (degrees) o f elements A , , and A , from E q . (5.12) 2 1 0 72 5.3c Inverse Fourier transform of the propagation matrix A The accuracy of the synthesized functions can be further checked by comparing their inverse Fourier counterparts in the time domain with those of the exact functions. Figure 5.5 shows the exact and synthesized functions of the first three elements of column 1 of a(t) (time-domain form of A(co)) for the line of Figure 4.2; the fitting accuracy is very high. The program of [55], with the Lanczos filter described in section 3.3a applied to the exact and synthesized functions, were used to obtain these inverses. Analogously to the suggestion in [36, p. 4-95], the physical meaning of a(t) can be explained (Figure 5.6) by connecting a set of currents I to node m, together with a shunt source characteristic admittance matrix Y (co), and shorting all phases at node k. In Figure 5.6, c summing the currents at node m gives I . Isource = Yc(co)V + m m x - 1 0 5 [ \ Exact function Synthesized function - L ^ \ \ a i i ( t ^--" 21 a Js^T^' 3i a 0.331 ) ( t ) ( t ) Time (msec) 0.34 Figure 5.5 - Time-domain propagation functions: a,,(t), -a (t), and-a (t) 21 31 73 Ik I k m m > source Shorted Figure 5.6 - Current source with matching admittance Yc(co) connected to node m Then, from Eq. (5.10), with the direction of the currents in node k the same as that from node m, Ik = A lsource If I source (5.22) is 1.0 for phase #1 and zero for the other phases, for all frequencies, then I bek comes the first column of the matrix A. This situation corresponds, in the time domain, to the injection into phase #1 of a unit current impulse at node m, with the other phases open. Therefore, the inverse Fourier transform of A,, must be the current impulse that arrives at node k. The inverses of the other elements in the first column of the matrix A are the currents of the other phases arriving at node k, induced through coupling between phases. It is seen in Figure 5.5 that the term a,,(t) has positive polarity and a relatively high magnitude, while a (t) and a (t) have negative and smaller magnitudes (-a (t) and -a (t) are plotted). 21 3! 21 31 These functions exhibit the typical attenuation and distortion of impulse waveshapes. 5.4 Synthesis of Y and A for asymmetrical transmission line configurations c Figure 5.7 shows a railway running parallel to a 1 mile long three-phase 138 kV line. The rails are represented as conductors and the ground wire is assumed to have uniformly zero potential. Thus, the line configuration is a five-phase system. This line configuration 74 6 ft G.W.: 1-159 MCM A C S R _4-^>, A D i a . = 0.576 in. 15 ft Phase conductor: 1-1033 MCM A C S R b Q<^ i - i u j j mum ftivSK D o < Dia. = 1.212 i n . 7 ftf ? 1 f Q C 14 ft t a Length = 1 mile 5 I #1 o Rails #1 & #2: D i a . = 9.19 in. #2 28.5 ft f t 40 ft I o< 10 ft V 100 Ohm-m ground resistivity Figure 5.7 - Railroad near a three-phase 138 kV line is extremely asymmetrical. Figures 5.8 and 5.9 show the elements of column 1 of its characteristic admittance and propagation matrices, respectively. The agreement between the fitted and exact functions is very good, similar to the case of the double-circuit line of Figure 4.2. Table 5.1 lists the average number of poles that the fitting routine FIT-S typically allocated for the synthesized functions. Table 5.1- Average number ofpoles allocatedfor synthesizedfunctions Function Diagonal elements Off-diagonal elements A 24 30 Y 14 17 . c 75 x 0 1 Elements order from top down: 1, 2, 3, 5 & 4 1 0 1 0 1 0 1 0 1 0 1 0 1 0 Frequency (Hz) 1 0 0 8 (D <U 0 6 01 La 0> w "5b c < a in 4 0 2 0 <a OH 0 2 0 1 0 1 0 1 0 Frequency (Hz) Figure 5.8 - Elements o f column 1 o f matrix Y c for an asymmetrical line 76 Exact function 0 . 8 L Synthesized function D i a g o n a 0.6^ -a 3 J2 0 .4 O f f - d i a g o n a l s 0 . 2 L 0 1 0 T Maximum error < 3.0 % 1 0 1 0 1 0 1 0 Frequency ( H z ) 1 0 0 5 0 u 60 (D w 0 C < a - s o -10 0 1 0 1 0 1 0 1 0 Frequency (Hz) Figure 5.9 - Elements o f column 1 o f matrix A for an asymmetrical line 1 0 77 5.5 Open/Short -circuit tests in the frequency domain As described in section 4.5, open- and short- circuit calculations are performed to check the accuracy of the synthesized functions. The open-circuit voltages and short-circuit currents are directly evaluated in phase co-ordinates here. Applying boundary conditions to Eq.'s (5.9) and (5.10), these voltages and currents at the receiving end m become V - «„ = Yc [A + A - ] " V , (5.23) I -^/ = 2[A-A- ]" Y V . (5.24). 1 m , , <¥ k 1 1 m c k If a voltage column matrix E is applied at the sending end of the line (node k), then s V = E ,or k s Vm-^^Yc'tA +A-'f'Es, (5.25) Im-^ = 2[A-A- ]" Y E . (5.26) 1 1 C S 5.5a Zero sequence excitation The zero-sequence source voltage phasors in this test are real and all have 1.0 p.u magnitude. Figures 5.10-5.13 show the responses obtained for both lines (Figures 4.2 and 5.7) in this test. The approximate and exact responses are in very good agreement. The responses of phases b and c (not plotted) for the line of figure 5.7 also agree very well. 78 5.5b. Unsymmetric excitation For the double-circuit line of figure 4.2, the unsymmetric source voltage phasors (E ) s are the same as in section 4.5b, namely E sl =E s4 = \Z0°, E =E =0, s2 s5 andE =E s3 = 1Z120 , 0 s6 whereas for the line of Figure 5.7, the following values were used: E sa =.\Z0°, E sb = lZ- 120°, E = 1Z120 , andE 0 sc s raim = Es raim = 1Z0 . Figures 5.14 - 5.17 show the responses obtained for this test. The agreement between the exact and synthesized functions is very good. 79 2 5 2 0. Exact function Synthesized function 1 5. 1 0. 5 . 1 0 ^ 1 0 "2 1 0 L 1 0 1 0 Frequency (Hz) 1 0 1 0 1 1 0 Frequency (Hz) 1 2 1 0 Frequency (Hz) Figure 5.10 - Receiving end open-circuit voltages (p.u. magnitude) for line o f Figure 4.2 80 0 . 6 . o3 w 0) in ro 0 . 4 0 . 2 . 1 0 Frequency (Hz) 0 . 6 °a 0 . 4 . CM C/) <D in ro 0 . 2 1 1 0 0 1 0 1 0 1 0 Frequency (Hz) 0 . 6 CD 08 0 . 4 CO in <D w ro 0 . 2 1 0 1 0 1 0 Frequency (Hz) 1 0 1 0 Exact function Synthesized function Figure 5.11 - Receiving end short-circuit currents (p.u. magnitude) for line o f Figure 4.2 2 5 Exact function 2 0 ro 1 Lx 1 0 cu to ro Synthesized function 5. - 2 1 0 Aim, 1 o 1 o 1 0 1 0 1 0 1 0 Frequency (Hz) 1 0 'ro Frequency (Hz) 1 0 CM 8 - 6 - 'ro a: 1 0 1 0 1 0 Frequency (Hz) Figure 5.12 - Receiving end open-circuit voltages (p.u. magnitude) for line o f Figure 82 1 2 a) ro 1 0 1 0 1 0 1 1 0 Frequency (Hz) 'ro 0C 2 0 1 5 . 1 0 1 0 1 0 1 0 0 Frequency (Hz) 2 0 1 5- CM 2 ro 1 o 1 0 1 0 1 0 Frequency (Hz) Exact function Synthesized function Figure 5.13 - Receiving end short-circuit currents (p.u. magnitude) for line o f Figure 5.7 83 2 5 2 0 1 5 1 0 Exact function Synthesized function CO 0) CO 03 5 - 1 0 - 2 1 0 1 0 1 0 1 0 Frequency (Hz) 1 0 LO o3 CM CO cu CO co 4 1 . 0 Frequency (Hz) 1 2 1 0. CO o3 oo CO CO CO CD Frequency (Hz) Figure 5.14 - Receiving end open-circuit voltages (p.u. magnitude) for line of Figure 4.2 84 0 . 8 Exact function 0 . 6 Synthesized function 08 co 0 4 CO CD 0 . 2 0 I 1 0 1 0 1 0 1 0 1 0 1 1 0 Frequency (Hz) 0 . 2 °3 CM CO CO CO CD 0 . 1L x: Q. 1 0 1 0 1 0 0 Frequency (Hz) 0 . 8 0 . 6t co o3 co CO CO CO CD 0 . 4L sz CL 0 . 2L 1 0 1 0 1 0 Frequency (Hz) Figure 5.15 - Receiving end short-circuit currents (p.u. magnitude) for line o f Figure 4.2 85 1 4 1 2 .1 0 ro Exact function Synthesized function 8 cu CO ro o. c J 1 o 1 0 1 0 1 0 1 0 Frequency (Hz) 1 0 1 0 . 1 0 1 0 1 0 1 0 1 0 Frequency (Hz) CN ro or 1 0 1 0 1 0 Frequency (Hz) Figure 5.16 - Receiving end open-circuit voltages (p.u. magnitude) for line o f Figure 5.7 86 1 2 Exact function 1 0 |- Synthesized function 8 CD 8! 6 l- CD ^ 4 2 I- 1 0 1 0 1 0 1 0 1 0 Frequency (Hz) 2 0 1 0 1 0 1 0 1 0 1 0 Frequency (Hz) 2 0 1 0 1 0 1 0 . 1 0 1 0 Frequency (Hz) Figure 5.17 - Receiving end short-circuit currents (p.u. magnitude) for line o f Figure 5.7 87 5.6 Time-domain simulation tests To demonstrate the accuracy of the proposed frequency-dependent line model, timedomain simulations for transient as well as steady-state conditions were performed. The results are presented in this section. For testing purposes, the equations of the developed frequency-dependent transmission line model were interfaced with the U B C E M T P version. MicroTran through subroutine CONNEC. This interface is described in Appendix C. Additional tests and field measurement comparisons for typical lightning events will be shown in Chapter 6. 5.4a Open-circuit test with zero-sequence energization In this test, all phases of the double-circuit line of Figure 4.2 are energized simultaneously with a 1.0 p.u. step voltage at time zero, with the receiving end of the line open. Figure 5.18 shows voltages at the receiving end of phase 1. Note that the response produced 0 5 1 0 Time (msec) 1 5 Figure 5.18 - Receiving end phase 1 open-circuit voltage (p.u.) 88 by the new line model agrees very well with those simulated with the frequency-domain model [41], and with the existing EMTP frequency-dependent line model [34]. It is worth emphasizing that the EMTP frequency-dependent line model using a constant transformation matrix at 5 kHz also gives a very reasonable response in this test, because the zero sequence modal transformation vector has very little, if any, frequency dependence. 5.4b Magnetically induced voltages from a three-phase 138 kV line Figure 5.19 shows a railway running parallel to a three-phase 138 kV line. As mentioned earlier, the rails were modelled as conductors. In order to investigate the possible magnetically induced voltages on the rails and to show the performance of the new line model,- a test condition is set up as illustrated in Figure 5.19. In this test, all three phases are energized simultaneously at time zero with a symmetric 60 Hz voltage source through a 190 Q source impedance. A l l phases, including the two rails, are shorted at the receiving end. Magnetically induced voltages will appear at the open end of the rails. Figure 5.20 shows the induced transient and steady-state voltages of rail #1. Table 5.2 lists the RMS values obtained with the compared models, as well as the exact solution obtained from a phasor solution. The results in Figure 5.20 and Table 5.2 show that the existing E M T P frequency- dependent model, with constant transformation matrices evaluated at two different frequencies, produces reasonable answers during the transient period, however, it leads to inaccurate results in the steady-state condition. The steady-state induced voltages produced by the new line model agree much closer with the exact solution. Note that even though the transformation matrix was calculated at 60 Hz, the existing frequency-dependent 89 G.W. t=0 Length = 1 mile 190 Ohms ^JPw)f\j)i38kv Rail #2 M Rail #1 i Figure 5.19 - Magnetically induced voltages on rails 1 - New model 2 - Existing model w. constant Ti at 5 kHz 3 - Existing model w. constant Ti at 60 Hz CU so > c CD '</> c CD 0 . 0 0 . 0 2 4 0 . 0 0 8 6 Time (msec) 1 0 1 5 2 0 2 5 Time (msec) 3 0 3 Figure 5.20 - Magnetically induced voltage on rail #1 5 4 0 90 line model, which uses a constant transformation matrix, still provided wrong steady-state solution at this frequency. The reason is because this line model used only the real part of the transformation matrix in its solution algorithm. This assumption is not required by the new line model. Table 5.2 - Steady-state induced voltages (RMS) Line Model Rail #1 Rail #2 Exact Solution 16.41 V 17.49 V New Model 16.27 V 18.19 V Frequency-dependent Model w. constant T, 20.75 V 23.03 V 11.32 V 15.98 V at 5 kHz Frequency-dependent Model w. constant T, at 60 Hz 91 CHAPTER 6 Simulation of Lightning Surges and Field Test Comparisons 6.1 Overview Unexplained system equipment failures often occur during heavy thunderstorms, for which power system engineers try to find out the causes with lightning surge studies. Most engineers use computer programs today such as the EMTP for such studies. Two new models have been developed for such programs, to represent transmission lines and transmission towers. This Chapter covers the application of these two models for the simulation of four typical lightning surge types, as outlined in Chapter 1. Sections 6.2 and 6.3 emphasize the improvements accomplished with these models, by comparing the simulation results with field tests and with field recordings. 6.2 Tower/ground wire surges When lightning strikes the tower or the ground wires, overvoltages are developed across the insulator strings. If these overvoltages are above the voltage which the insulators can withstand, flashovers to the phase conductors will occur. Simulations of these scenarios have indicated that the line insulator voltages are influenced by the models used in representing the transmission towers. In the past, the transmission tower has often been repre- sented as a transmission line with uniform surge impedance. As pointed out in section 2.3, 92 the transmission tower surge impedance does vary as the wave travels from the top to the base. Thus, a nonuniform transmission line model is required. The simulations of the measured insulator voltages from the published experiments of Ishii et al [18] are analyzed in this section. double-circuit towers. They performed measurements on 500 kV Their experimental set-up is similar to the case where lightning strikes the transmission towers, except that the current source was low in magnitude. The circuit parameters for EMTP simulations are displayed in Figure 6.1. The circuit consists of three sections of a double-circuit line with the geometry of Figure 6.2, of three transmission towers, and of a low current source with a 400 Q source impedance. A 300 m cable with a surge impedance of 50 Q is connected between the source and the tower top. All of the conductors are terminated with 350 Q resistors at both ends. The rest of the parameters used are similar to those described in [18], with the exception of the line model: the authors of [18] used the Semlyen frequency-dependent line model, while a constant-parameter line model is used here. This difference has no noticeable effect on the insulator voltages in this case. The middle tower No. 7 is modelled as 4 transmission line sections (Figure 6.3), three uniform and one exponential -line, with the crossarms ignored. In concurrence with the measured impedance values in [18], the uniform sections have a uniform 220 Q surge impedance, while the exponential one has a nonuniform surge impedance of 220 Q at the upper end and 150 Q at the bottom end. Capacitances of the insulator strings and stray capacitances from the phase conductors to the tower are also modelled. The simulated and measured waveforms of the three crossarm insulator voltages (potential difference between crossarm and phase conductor) at tower No. 7 are shown in Figure 6.4. A time delay of 93 approximately 1 ps is seen in these voltages because of the connected cable between the current source and the tower top. The voltages are very close compared to the measurements. Thus, the simulation results supports the validity of representing transmission towers with the exponential transmission line model! It should be mentioned that the capacitances (from the insulator strings and from the phase conductor to the tower) influence and retard the electromagnetic field propagation in the transverse direction. This aspect is not taken into account by the conventional timedomain transmission line models because they only represent the T E M mode in the horizontal direction [56]. Without including these capacitances, the lower-phase insulator voltage would have the highest magnitude. This would be the opposite of what had been measured. No. 8 449 m 400m Constant parameter line Constant parameter line Uniform line model < . . ' Uniform line model 17 Q Exponential line model 17 n Figure 6.1 - Circuit for E M T P simulations 94 ^ 13-2 m (No. 56) 11.1 m(No. 6-7) 133m(No.7«) ^ Ground wire, 1.75 cm diameter If 19.0 m y 1 127m A 20.8 m o 2 O 127m 3 19.0 m 20.8 m O 4 Phase conductor 4 X 3.84 cm diameter at 50 cm spacing -5* O 5 o—o 6 O 6 ! o o 40X)m(No.56) 26X)m(No.6-7) 36.0 m (No. 7-8) ^ i 50 cm 500 Ohm-m ground resistivity Figure 6.2 - Conductor geometry at average height L i n e #1 Uniform Line #2 Uniform Line #3 Uniform CR#1 CR#2 1 CR#3 Line # 4 Nonuniform (Exponential line) / / / / / / / / / / / ~~" "~ " Figure 6.3 - Tower representation 95 Measurement Simulation 1 2 0 40 V / d i v . 50 « (89%) Time (us) 1 2 0 > 0 1 2 3 Time (us) Figure 6.4 - Simulation of Ishii et al experiment and measurement from [18] 4 5 96 6.3 Simulation of a recorded shielding failure event in a 765 kV system A simplified diagram of the American Electric Power (AEP) 765 kV system surrounding the Rockport area is shown in Figure 6.5. A monitoring system has been installed at Rockport since 1986 [2]. This system has captured many interesting shielding failure events during heavy thunderstorms. Shown in Figure 6.7 are the recorded waveforms for one such event. The expanded waveforms around the strike location are shown in Figure 6.8. A careful examination of these waveforms revealed that a shielding failure had occurred on phase A of the Rockport-Sullivan line, approximately 31 miles away from Rockport. The travelling surge reflections denoted by a-d in Figure 6.8 are identified on the Rockport phase A voltage as: • (a) - Initial lightning surge arrives at Rockport, which also produces coupled surges on phases B and C. • (b) - First reflection from Sullivan comes. The travel time between this and the first surge of (a) suggests that the stroke hit the Rockport-Sullivan line, approximately 31 miles from Rockport. • (c) - Second reflection from Sullivan arrives. • (d) - First reflection from Jefferson arrives. In the simulation of this shielding failure event, all of the transformers are ignored, except for the one at Jefferson station, The transmission lines, system linear reactors, and system 60 Hz sources behind their appropriate short circuit impedances are represented. The struck line (Rockport-Sullivan) was first represented with the new transmission line model and then with the existing frequency-dependent model. Because of insignificant 97 contributions to the simulation results, all other secondary lines were modelled with untransposed constant-parameter lines, calculated at 5 kHz. With the exception of the Rockport-Jefferson line, these secondary lines do not have to be included. The reason for including them in this simulation is because the author wanted to investigate if there are any significant reflections coming back to Rockport from far away stations. It was concluded from observing the results that no significant reflections arrived back to Rockport from these far away stations. According to calculations performed in Appendix D, the maximum current for lightning strokes which can escape the shielding system (the two ground wires) and terminate directly onto the phase conductors is 8.13 kA. The event was simulated by injecting a current source with this magnitude into phase A at the strike location, on the Rockport-Sullivan line. The injected current waveshape had the concave wavefront recommended in [7]. Preliminary results showed that the struck phase voltage was noticeable lower in magnitude than that of the recording. Subsequently, the current source was increased to 11.2 kA. The other characteristics of the current source are listed in Figure 6.6. The equations for the current source are: i(t) = -{At (t/t ) + Bt (t/t ) } n 90 90 90 90 ]e (front), f (tail), 90 /( ) = -{/ -('-'9o)/'._/ g-(W9o)/'2) ? for t<t 2 or >t t 90 where the front part has At^ = 5324.225, BV = 4755.775, t^ = 31.034483 [is and n = 8.0163697, and the tail part has I, = 11571.533,1 = 1491.5327, t, = 112.98168 \xs and t = 2 2 0.99275381 |xs. With this current source, the simulated waveforms are shown in Figure 6.9. As can be seen from these plots, the waveforms are in very good agreement with the recordings of 98 Figure 6.8, with respect to the magnitudes, travelling times and polarities. It can also be observed that the initial surges coupled into phases B and C produced by the new transmission line model match the recordings closer than those given by the constant T, frequency-dependent line model. It should be noted that a current magnitude of 11.2 kA is as probable as the value of 8.13 kA calculated with the electrogeometric theory in Appendix D. This electrogeometric model is only an approximation, and has a margin of error. In addition, it has been reported that subsequent strokes in a single flash sometimes have higher magnitudes than the first stroke (with a factor of up to 200%) [44]. The event described in Figure 1.5 is evidence of this, where the third stroke actually caused a flashover, which most likely had a higher magnitude than the earlier strokes. 99 36.5' 36.5' 18"^ A h c A 34.75" K— 6 6 o o o o L A .W. 5^" 44.5' X 4 4 5 ' > hq = 99.5' =57.5' Line geometry at average height D U M O N T S T R I K E L O C A T I O N AAAA H A N G I N G R O C K R O C K P O R T Figure 6.5 - Simplified diagram of 765 kV system surrounding Rockport A rest "11-2 kA Tfront = 20 LIS T =100Lis S = -1.4 kA/us C tail max \ - 6 0 0 0 - 8 0 0 0 1 0 0 0 0 1 2 0 0 0 Note that T , is defined throughthe 30% and 90% points [7]. V \ 2 0 flon 4 0 Time (us) 6 0 Figure 6.6 - Lightning current source (A) 8 0 1 0 100 x o O) s? 0 1 0 . '5 U o > •e oQ. 0 . 5 J< O O 1 o > k 0 . 5 U m © ra o a 0 . 5k o o a: 0) o> 0 1 V . 5 L o > (J a> in re t o a 0 : 5 O o a. Figure 6 . 7 - Lightning shielding failure waveforms recorded on Oct. 7 , 1990 102 103 6.4 Shielding-failure flashover If currents in the order of tens of kA or less escape the shielding and terminate directly onto the phase conductors, it is most unlikely that these shielding failures will flash over. This is not always the case, however. Shielding failures can result in flashovers, provided that there is a combination of worst scenarios. For example, considering the 765 kV system of Figure 6.5 again, if a subsequent stroke of approximately -16 kA in magnitude is injected at a negative peak of the system steady-state voltage, a voltage larger than the line standard BIL would appear across the line insulators. For a normal 4.4 m 765 kV line insulator, a value of 2600 kV critical flashover voltage is assumed for lightning strokes having a 20 ps wave front [5]. The lightning current source in this case has a similar waveshape with the one in Figure 6.6, with A crest = -16 kA, T front = 20 ps, T tail = 100 ps, S max = -2.0 kA/ps. A vol- tage-dependent switch is used between phase A and ground to model the flashover. This switch closes as soon as the voltage across it is greater or equal to the 2600 kV BIL. Figure 6.10 shows the simulated phase-to-ground voltages at Rockport when the current source was applied 2 miles away, on the Rockport-Sullivan line. Note that Phase A shows a voltage close to zero during the fault period (about half a cycle). The purpose of the foregoing analysis is to identify situations where shielding-failure flashovers could occur. It is not within the scope of the present thesis to investigate the details of the breakdown mechanism of the line insulation. There is a large volume of experimental work on this topic [5, 7, 44]. To assume that flashover occurs just when the voltage exceeds the breakdown volt-time characteristic is an approximation, which is useful for preliminary studies. More detailed line insulation co-ordination studies can then be done as outlined in 104 the above references. Experimental volt-time characteristics are only valid for standard laboratory impulses. For non-standard lightning waveshapes, as typically seen in field recordings, these characteristics are of limited use in predicting the performance of lines. 6.6 Back flashover The circuit of Figure 6.11 is used for simulating back flashover events. It represents a 765 kV line of the geometry shown in Figure 6.5, and consists of 13 towers with each span being 0.25 miles. To avoid reflections coming back, the line is terminated with its own surge impedance at both ends. The untransposed constant-parameter line model, calculated at 500 kHz, is used for the line spans between the towers. The struck tower (tower #7) is modelled as an exponential transmission line, while the others are modelled as constant-parameter lines. With 60 Hz voltages connected to both ends of the line, the lightning current source shown in Figure 6.12 is injected at the tower top of tower #7. Since phase A was assumed to be at its maximum 60 Hz positive voltage at time zero, its line insulator will see a higher voltage compared to other two phases. Thus, a back flashover will most likely occur to this phase first. Figure 6.13 shows the tower top, phase A , and phase A insulator voltages obtained at tower #7, for this incident. Note that the phase A insulator voltage reaches a maximum value of -1.7 M V . With the volt-time characteristics of [5], for a 4.4 m 765 kV insulator, a critical flashover voltage of 3.1 M V is obtained for standard lightning waveforms having a 3 us wavefront. Therefore, for a flashover to develop, the injected current must be increased by a factor of (3.1/1.7). Thus, a back flashover could occur if a current of 45(3.1/1.7) or 82.6 kA magnitude, with the same front time, hits tower #7. The same back 105 flashover scenario can be established by increasing the tower footing resistance to a larger value, without having to change the magnitude of the lightning current. Figure 6.14 shows the respective voltages when the tower #7 footing resistance is increased to 180 Q. Thus, one can conclude that back flashovers could occur on E H V transmission systems provided that the lightning current magnitude is large enough, or that the tower being hit has a high footing impedance. It must be noted, however, that E H V towers are usually designed with low footing impedances, with the exception of those in rocky mountainous locations. For transmission lines over this type of terrain, simulations are important to predict the line outage rates, and to analyze the influence of transients on protection systems. Summary Lightning studies for tower/ground wire surges and back flashover events are only concerned with a few tower spans close to the strike location, because the surges attenuate quite fast in the ground (zero sequence) mode as the distance increases. For shielding-failure and shielding-failure flashover events, however, the entire struck line has to be modelled because the waveforms may still be significant at a far distance down the line. For example, referring to Figure 6.5, it has been observed that a reflection coming back from Sullivan reinforced the oscillations which were initially created by the shielding failure surge; in the high-voltage transformer tertiary windings at Rockport [57]. Thus, even after travelling the length of the line, a reflected surge could still excite the transformer natural frequency oscillations, and result in overvoltages in the windings. For this reason, it is best to represent all incoming lines to the substation with their entire length, in investigating the effects of 106 shielding failure and shielding-failure flashover waveforms impinging on substation equipment. For most practical purposes, the struck tower can be represented with a one-section exponential line model, as shown in Figure 6.11, which gives good enough results and is easy to use in simulations. For the situation where the insulator voltages of different crossarms need to be known more accurately, the transmission tower must be broken up into and modelled as a number of exponential line sections. The number of sections depends on the number of crossarms which the tower has. For the situations of tower/ground wire surges and back flashover, the line representing each individual span between towers should include the ground wires as additional conductors, and represent the towers along the path, whereas in shielding-failure and shielding-failure flashover cases, the line can be modelled with the ground wires eliminated. The accuracy of a particular line model affects greatly the results in those events which involve the phase conductors directly, such as the shielding failures and the shielding-failure flashovers, whereas it has less influence in the tower/ground wire surges and back flashover cases. The reverse holds true for the transmission tower representation. 107 108 Exponential line Z = 300 and 100 n at top and bottom, respectively. c # U n i f o r m line Z =300fi c Figure 6.11 - E M T P circuit for simulation o f a back flashover event X 10 A crest = -45.0kA Tfront = 3 U.S T,ail = 75 JXS S max = -37.5 kA/ns Figure 6.12 - Lightning current source ( A ) 109 Figure 6.13 - Simulated voltages at tower #7 with a 15 Q footing resistance 110 x 0 1 0 6 . 5 I Time (us) Figure 6.14 - Simulated voltages at tower #7 with a 180 Q footing resistance Ill CHAPTER 7 Future Research in Transmission Line and Transmission Tower Modelling Although the main goals set for this thesis have been accomplished, the author believes that improvements can be made if the following additional research is carried out. 1) Extend the direct phase-domain modelling technique of overhead transmission lines to underground cable modelling. The method described in this thesis is general and could potentially be applicable for the case of underground cables. 2) Include corona effects in the simulation of lightning surges, if reliable corona models are available. This aspect has been ignored in this thesis due to the lack of field test verifications. It is appropriate to note here that the upper part of the wavefront is most likely being reduced by corona, as the surges travel along the line. This wavefront reduction has been noticed in the recorded waveforms in reference [2]. This system is in a remote location, and has been experiencing frequent lightning incidents 50 miles or farther away. The recorded wavefronts have been somewhat slow (> 20 ps), compared to typical lightning data in the literature. The corona effect is thought to have a major influence on this. 3) Develop equations for the calculation of the surge impedances of the exponential line model, which can be used to represent typical high voltage transmission towers. The tower surge impedance is a function of its dimensions. When the conductor geometries are specified in line constants programs, the tower dimensions are indirectly given as weil, and it may 112 therefore be possible to modify line constants programs to produce tower surge impedances as well. Built-in models in such programs should cover most of the typical towers presently used in industry. The common classes are described in [5, 15], which includes the lattice, lattice and steel pole, H-frame, guyed-Y and typical 735-800 kV towers. 4) Modify the E M T P to reduce storage requirement for cases involving long transmission lines and A C steady-state initial conditions. The derivation of this modification is shown in Appendix E . 113 CHAPTER 8 Conclusions E M T P simulations have become widely accepted for lightning transient studies, which are done to prevent equipment failures and transmission line outages. In such studies, it is important to realistically simulate and understand all possible events created by lightning. From preliminary EMTP simulations for these events, it was seen that the representation of the transmission tower and transmission line influences greatly the simulation results. It was also observed that the existing models used in representing towers and lines have limitations. For this reason, new and improved models were developed in this thesis. A n exponential tapered line model was derived for the tower. Its time-domain simulation results showed good agreement with those produced by a cascaded model where the line is divided into many uniform sections. The model is numerically stable because it syn- thesizes the line functions in the frequency domain with stable rational functions of minimum-phase-shift type. It includes these functions in the time domain by means of recursive convolutions. The model is suitable for the representation of transmission towers because the tower, from a travelling wave standpoint, is a nonuniform structure. It can also be used to model general nonuniform transmission lines with complicated space-dependent characteristic impedances, as described in [21]. A simple, fast and accurate method based on the Newton-Raphson iteration technique was implemented for finding transmission line modal parameters. It overcomes the . • • . • • 1 1 4 eigenvector and eigenvalue switching problems experienced by most other methods. It produces modal parameters which can reasonably be fitted with rational functions of minimumphase type. This should improve the modal transmission line formulation of [35]. A new wideband time-domain transmission line model was developed. It includes the complete frequency-dependent nature of untransposed overhead transmission lines directly in phase co-ordinates. It. avoids the assumption of a constant transformation matrix, which the existing frequency-dependent transmission line model uses. Extensive simulation results were presented to show that the new model provides reasonably accurate answers for both steady-state and transient studies. The model is numerically stable and can be used for most types of transient studies. In terms of simulation time, the proposed model was, on average, about 25% slower than the existing frequency-dependent model using a constant T,. However, this could be improved if the programming code is optimized. Typical lightning surges were simulated to compare the solutions produced by the developed models with actual measurements. With these new modelling approaches, simulations matched field recorded transient waveforms closer than existing models. 115 APPENDIX A Asymptotic Behaviour of Modal Transformation Matrices at Very Low and High Frequencies It is shown in this appendix that the elements of the matrix T, of Eq. (4.3) become real and approach the asymptotic values at the high and low frequency limits. A . l Behaviour at very low frequencies When the frequency approaches zero, the conductor internal impedance matrix Z bec comes the d.c. resistance matrix Rj,.. Then YZ=;coP-'[/-coLi s P + Z + R c], 0 0 e (Al) d ignoring the shunt conductances of the line. Since the first two terms in the right-hand side of the above equation approach zero faster at very low frequencies, it follows that Y Z —> ycoP Rd as co - » 0. With P" is the _1 1 C capacitance matrix C , the Y Z product of Eq. (Al), at low frequencies, is YZ«/coCR . d c (A2) Because the eigenvectors stay unchanged if the matrix to be diagonalized is scaled with a constant, Eq. (A2) can be reduced to just the product of CRdc- If all of the conductors have the same dc resistance value, the matrix to be diagonalized is then reduced further to just the matrix C . Thus, the eigenvectors of the Y Z product at very low frequencies can be found by diagonalizing only the matrix C , and Since C is real and symmetric, its eigenvectors are real and constant. This agrees with the solution method described in [33]. 116 A.2 Behaviour at very high frequencies At very high frequencies, the elements of the two matrices Z a n d Z in Eq. (4.6) become c pi j(o\ifd i (A3) e 2n V r )> 2nrd,el ycop 2n 2d (yi +yj) e2 0 E (A4) Dl Eq. (A4) is obtained because \d \ << Dy, for co -> oo, and e2 In WA ^D] + Ad iy +yj) j e2 i «ln Ad {yi +yf) 1+ e2 = ln Dl 2d (yi +yj) e2 because for small 0, ln(l + 0) —» 0. 4 Then, Y Z = -w u s U + 2 0 0 1 p - . B 2TIZ (A5) 0 where U is the unit matrix, and the elements of the matrix B are: B =^ ii + \i d ^andB r 2d (yi +yj) e2 e i J D (A6) where u. = 77- is the conductor relative permeability. r Because d \ and d e e2 -> 0, as (a ->• 00, the elements of the matrix B of Eq. (A6) event- ually become negligible at very high frequencies. Thus, from Eq. (A5), Y Z - » - c o p s U , asm -> 00 , 2 0 0 (A7) Therefore, the eigenvectors of the Y Z product approach constant and real values as the frequency increases to infinity. 117 Eq. (A7) is also obtained in [36, p. 4-58] by making the two assumptions: 1. Conductor resistances and ground return resistances are ignored. 2. At very high frequencies, all currents flow on the surface of the conductors, and on the surface of the earth. This is because the complex depth of penetration constants d el and d approach zero at high frequencies. e2 In actual calculations, losses are small, and that the complex depth of penetration constants are not exactly zero at the 1 MHz limit. Therefore, the eigenvectors have very small angles (close to zero) at this limit. This has been seen in all the eigenvectors obtained for the lines tested in this thesis. 118 APPENDIX B Recursive Convolutions B l . In matrix form If the elements of Y are approximated by rational functions with real negative poles and c zeros, then the numerical convolution of the y (t) matrix and the voltages vector v(t) has the c following form, y c ( 0 * v ( 0 = g(0 = y e q v ( 0 + M O , where is the history current vector evaluated from g(t- At) and \(t- (Bi) At), and y e q is a real, constant symmetric matrix. In a similar fashion, the elements of the A(co) matrix can also be synthesized with rational functions, but with a time delay associated with each element. Numerical convolution of the corresponding time-domain function a(t) with a forcing vector f(t) is evaluated as a(0 * f(0 = r(7) = W ( 0 The currents vector \ p(f) depends exclusively pr0 (B2) on past values of f(t) and r(t). Introduc- ing (Bl) and (B2) into (5.13) and (5.14) leads to (5.17) and (5.18), with = iprop(t) - hi (0 • (B3) B2. Integration constants and scalar convolutions When rational function syntheses are performed, each element of the matrix y (t) bec comes a simple sum of exponential functions and contains a constant, namely 119 m yc-ij(t) = k d(t) +1 kie-P*u(t). 0 (B4) i=\ In a slightly different form, each element of the a(t) matrix contains a unique and different time delay Xy associated with it, and the constant k becomes zero in this case, namely 0 m = k,e-P>^ (B5) i=\ Consider the matrix convolution of Eq. (Bl). This convolution operation is carried out as a normal matrix multiplication, but each multiplication involves a scalar convolution [30, 34, 35]. For example, the Convolution between element 11 of the y (t) matrix and element 1 c • of the v(t) matrix can be performed recursively as y -i i (0 * v, (0 = (t) = dv , (0 + hi(t), c g (B6) where m d = k + zZdi, 0 (B7) /=i m h (t) = T,c g (t-At) + e vi(t-At), i-\ i i i i (B8) , c di and e; are integration constants which depend on the integration method. For the trapi5 ezoidal rule of integration, they are defined as c, = ! ^ £ ; .*=., = - ^ £ - . 2+piAt (B9) 2+pjA It should be noted that the stability of the trapezoidal rule of integration was proven and has been used in most of the present EMTP models [36, 34, 35]. The scalar convolution of an element of the matrix a(t) and a forcing function from equation (B2) is also evaluated recursively, but based solely on past values, m a j(t)*fx(i) = r(t) = zZr (i), i i (BIO), 120 where n(f) = cir,(t- At) + dfi (t- Tij) + e(f\{t- Ty - At). The integration constants for Eq. (Bl 1) are the same as those from Eq. (B9). (B11) 121 APPENDIX C Interfacing User-defined Models with the EMTP version MicroTran The equations of the developed models in this thesis have been interfaced with the E M T P version MicroTran, in the time domain, through the subroutine C O N N E C . This is a convenient way to test new models. The final implementation will be done directly in the main time step loop of the program, in a much more efficient manner. Consider a transmission line connecting node k to node m in Figure CI. MicroTran can describe the behaviour of the network seen by both terminals of the line, in the form of a Thevenin equivalent circuit, with the equation Vk Vko R THEV-k 0 hm v Vmo 0 RTHEV-M Imk m (CI) The Thevenin equivalent resistance matrix has zero off-diagonal elements if the rest of the network has no lumped-element connections between nodes k and m, which was true in all cases studied in this thesis. k m ex km 'mk THEV-k m ^THEV-m 9 m Figure C I - Equivalent circuit seen by the transmission line terminals 122 At each time step, the subroutine C O N N E C supplies values for the open-circuit voltages and the Thevenin equivalent resistance network. These parameters, together with the equations of the line are used to solve for the line terminals' voltages and currents. The equations are expressed as follows: C.l Single-phase Exponential transmission line equations From Eq.'s (3.21) and (3.11) derived in section 3.2b of Chapter 3, ikm Zkeq 0 imk 0 Z ep ^mep -1 Vk Vm m m v eh-k (C2) @ h-m <^h-m Substituting Eq. (C2) into Eq. (CI), . v k v m R THEV-k Vko 0 V mo -1 0 Zkeq RmEV-m 0 0 Vk - Zmep Vm eh-k eh—m Hence, 1 0 0 1 RmEV-k 0 0 Zkeq 0 0 zmep RmEV-m Vko Vmo + R THEV-k 0 0 RmEV-m Vk V m Zkeq 0 eh-k 0 z,mep e/i-m (C3) The only unknown variables in Eq. (C3) above are the voltages v and v . Thus, these k m voltages can be solved and obtained at each time step. The currents can then be found from Eq. (C2). C.2 Multiphase frequency-dependent transmission line equations For a multiphase transmission line system as in Figure 4.1, Eq. (CI) contains matrix quantities instead, 123 v* v 0 RTHEV-IC 0 0 R'77/£K-»i m (C4) From Eq.'s (5.17) and (5.18) derived in section 5.2 of Chapter 5, yeq . V w 0 - -l 0 y 1km + lhist-k (C5) Imk + ^hist-m e ? Substituting Eq. (C5) into Eq. (C4), yeq Urn+ 0 0 y e q y.eq 0 0 Yeq imk lhist-k Vko R-THEV-k lfiist-m Vmo 0 0 ^THEV-m Hm imk Hence, J yeq 0 0 -1 _L RmEV-k yeq _ Thus, the currents i found from Eq. (C5). 0 k m Vyt 0 Rj-ffiF-m and i m k imk -1 0 0 yq _ e (C6) lhist-m are solved first with Eq. (C6), the voltages can then be 124 APPENDIX D Estimation of Maximum Lightning Shielding Failure Stroke on a Typical 765 kV Line Based on Electrogeometric Theory G.W. o\o Conductor Bundle ~5~oo A h = 129.5' - (40')2/3 = 99.5' = 30.33 m g h c = 99.5'- (50')2/3 = 57.5'= 17.53 m a = 2.44 m d=13.03 m a=l0.£ > Figure D l - Shielding angle of the right conductor bundle Figure D l shows the location of a ground wire and the closest phase conductor bundle at their average heights for a typical three-phase 765 kV flat transmission line. The average height is equal to the height at the tower, minus two thirds of the sag. In this conductor arrangement, shielding failures can only occur with the two outer phases. According to reference [5], the following formulas are used, s=\or°- , 65 [3 = 0.8, where S is the striking distance and [3 is the correction factor. The meanings of these parameters are illustrated in Figure D2. It is observed that, based on the striking distance theory, 125 the conductor bundle is vulnerable when lightning enters the uncovered area X . The striks ing distance S is determined by the lightning stroke magnitude I. Thus, if h , h and a are g c held fixed, the uncovered distance X would be diminished when the current I increases until s reaching its maximum value. With this condition, Figure D2 is modified to D3. From Figure D3, PSmax = 0.8S , max = h +S , c -> sinG = 0.8 max sin9 , 17.53 (Dl) Also sin 0 = sin 1 0 0 . 8 ° - cos" (^-] 1 Expanding by using the trigonometric identity, sinO = sin 100.8° • cos _ / 6.5 ^ v -> sine 6.38 - cos 100.8° •sin \s~J c o s "J m a x v "->max ' y -cos 1 0 0 . 8 ° - s i n cos _ , f 6.5 "l — ; v >->max (D2) ' Combining Eq.'s (Dl) and (D2), 6 7^-cosl00.8°sin <J cos _i( 6.5 ^i->max max = 0.8- 17.53 ' Rearranging, cos 1 0 0 . 8 ° s i n i c o s 23.91 _ / 6.5 0.8 = 0. Dividing Eq. (D3) by coslOO.8' —> sin cos if 6.5 } + ^ ^ - 4 . 2 7 = 0. Solving for S„ S max = 38.83w. (D3) 126 Since S m a x =10»/° Therefore, 7 max -> 6 5 / = 0.029 • = 0.029(38.85) 154 S , lM = 8.13 It is observed from the above procedure that the maximum lightning current which can penetrate this particular line shielding is 8.13 kA. This value is comparable to that obtained from the modified equations suggested in [44], which were also derived from the electrogeometric model of [5]. The concept of attractive radius was included in these modified equations. For helpful purposes, they are carried out here. The equation used for the current is still the same as that listed above, which is 7 = 0.029^ , 5 4 max but with S calculated as f -B - JB 2 S S=Y +A C S S a V where: =^ n =23.93, Xc-x 2.44 g = h^r m , n 30.33- 17.53 = 0 - 1 o 9 n A 0 6 ' p = 0.36 +0.168 l n ( 4 3 - F ) = 0.85 , 0 A = m -m p -$ = 2 2 2 2 s -0.7213, C = (m + \)= 1.036, 2 S B = p C , = 0.886. S Hence 5=35.84 -> /max = 0.029(35.84) 154 = 7.18 kA. Striking distance Figure D2 - Incomplete Shielding, Width X is uncovered s Maximum striking distance GW Conductor Bundle ps, max a = 79.2° 6 = 180°-a-co = 100.8° - cos- (z/2S ) 1 max SlL. Figure D3 - Effective shielding, Unprotected width X is reduced to zero s 128 APPENDIX E Representation of a Long Line with Linear AC Steady-State Initial Conditions When the maximum simulation time t ^ is less than the travel time x of a long line, the line will' only be represented by its surge impedance matrix at the two ends if the case starts from zero initial conditions in the MicroTran version of the E M T P [36]. Computer time and storage space are then significantly saved because there is no need for recording the history terms for such a long line. This simplified model is most often used in lightning studies. However, one often wants to include the system steady-state voltages in lightning simulations as well. In that case, the simulation has to start from linear A C steady-state conditions. This can be accommodated by a simple modification, by adding A C steady-state history current sources at both ends of the line. This simplified modification can be illustrated as follows: For simplicity, consider a single-phase transmission line in Figure E l . Its current and voltage relationships at node k can be expressed as v*(0 + histk, ikm(J) = km (El) Z mk ) Figure E l - Single transmission line m 129 where hist = k V m ( f ~ i ) -i (t-x). (E2) mk Similarly, the same expressions can be applied for node m, but with k and m subscripts interchanged. With zero initial conditions, the line can then be replaced by shunt resistances with a value of Z at both ends, as shown in Figure E2. Figure E 2 - Equivalent circuit o f a long line without initial conditions For linear A C steady-state initial conditions, additional current sources are needed, thus this equivalent circuit is modified to the circuit of Figure E3, where hist is the A C current k source of Eq. (E2). The same expression with subscripts k and m interchanged is applied for hist™. m Figure E3 - Equivalent circuit o f a long line with linear A C initial conditions 130 REFERENCES [I] H.W. Dommel, "Digital Computer Solution of Electromagnetic Transients in Singleand Multiphase Networks," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-88, No. 4, pp. 388-399, 1969. [2] J.M. Schneider, E.N. Fromholtz, D . K Nichols, B.J. Ware, "The Rockport Transient Voltage Monitoring System," CIGRE Report 23-04, 1988 Session. [3] A. Greenwood, Electrical Transients in Power Systems, 2nd edn., Chapter 10, Wiley (1991). [4] K. Ragaller, Current Interruption in High-Voltage Network, Plenum Press, New York (1978). [5] J.D. Anderson, "Lightning Performance of Transmission Lines," Chapter 12 of Transmission Line Reference Book-345 kV and Above, 2nd edn., EPRI (1982). [6] CIGRE WG. 33.02, "Guidelines for Representation of Network Elements When Calculating Transients," CIGRE Technical Brochure, 1990. [7] CIGRE W G 33.01, " Guide to Procedures for Estimating The Lightning Performance of Transmission Lines," CIGRE Technical Brochure 63, 1991. [8] P.P. Barker, et al, "Characteristics of Lighting Surges Measured at Metal Oxide Distribution Arresters," IEEE Trans, on Power Delivery, Vol. 8, No. 1, pp. 301-309, 1993. [9] M.D. Perkins, et al, "Summary of Transient Data Obtained from Long-Term Monitoring of a 138-kV and a 500 kV Transmission Line," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-103, No. 8, pp. 2290-2298, 1984. [10] K. Nakamura, et al, "Artificially Triggered Lightning Experiments to an E H V Transmission Line," IEEE Trans, on Power Delivery, Vol. 6, No. 3, pp. 1311-1318, 1991. [II] C F . Wagner and A.R. Hileman, "A new Approach to The Calculation of The Lightning Performance of Transmission Lines III-A Simplified Method: Stroke to The Tower," AIEE Trans, on Power Apparatus and Systems, Vol. 79, pp. 589-603, 1960. [12] M . Darveniza, et al, "Modelling for Lightning Performance calculations," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-98, No. 6, pp. 1900-1908, 1979. [13] C F . Wagner and A.R. Hileman, "A new Approach to The Calculation of The Lightning Performance of Transmission Lines III-A Simplified Method: Stroke to The Tower," AIEE Trans, on Power Apparatus and Systems, Vol. 79, pp. 589-603, 1960. 131 [14] M . Kawai, "Studies of The Surge Response of a Transmission Line Tower," IEEE Trans, on Power Apparatus and Systems, Vol. 83, pp.30-34, 1964. [15] P.G Buchan and W.A. Chisholm, "Surge Impedance of H V and E H V Transmission Towers," Final Report for Canadian Electrical Association Contract #78-71, 1980. [16] M . A . Sargent and M . Darveniza, "Tower Surge Impedance," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-88, No. 5, pp. 680-687, 1969. [17] W.A. Chisholm et al, "Lightning Surge Response of Transmission Towers," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-102, No. 9, pp. 3232-3242, 1983. [18] M . Ishii et al, "Multistory Transmission Tower for Lightning Surge Analysis," IEEE Trans, on Power Delivery, Vol. 6, No. 3, pp. 1327-1335, 1991. [19] L.O. Barthold, "An Approximate Transient Solution of the Tapered Transmission Line," AIEE Trans., PAS 77, pp. 1556-1561, 1958. [20] C. Menemenlis and Z.T. Chun, "Wave Propagation on Nonuniform Lines," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-101, No. 4, pp. 833-839, 1982. [21] M . M . Saied, A.S. Alfuhaid, M . E . ElShandwily, "S-Domain Analysis of Electromagnetic Transients on Nonuniform Lines," IEEE Trans, on Power Delivery, Vol. 5, No. 4, pp. 2072-2083, November 1990. [22] H.V. Nguyen, H.W. Dommel, J.R. Marti, "Tower Models for Lightning Surge Simulations," IEEE Paper No. 045-5 PWRD, Winter Power Meeting, New York, New York, Feb. 1994. [23] E . A . Oufi, A.S. Alfuhaid, M . M . Saied, " Transient Analysis of Lossless Single-Phase Nonuniform Transmission Lines," IEEE Trans, on Power Delivery, Vol. 9, No. 3, pp. 1694-1701, July 1994. [24] M.T. Correia de Barros, M . E . Almeida, "Computation of Electromagnetic Transients on Nonuniform Transmission Lines," IEEE Paper No. 397-0 PWRD, Summer Power Meeting, Portland, Oregon, July 1995. [25] W.A. Chisholm and Y . L . Chow, "Travel Time of Transmission Towers," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-104, No. 10, pp. 2922-2928, 1985. [26] Working Group on Lightning Performance of Transmission Lines, "A Simplified Method for Estimating Lightning Performance of Transmission Lines," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-104, No. 4, pp. 919-927,1985. [27] J.T. Whitehead, "Lightning Performance of TVA's 161-kV and 500-kV Transmission Lines," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-102, No. 3, pp. 752-768, 1983 132 [28] W . A ChisKolm, "Lightning Surge Response of Transmission Towers and Lines," Ph.D. Thesis, University of Waterloo, Department of Electrical Engineering, 1983. [29] R.S. Elliott, A n Introduction to Guided Waves and Microwave Circuits. Prentice-Hall, Inc., [13] M.A. Sargent and M . Darveniza, "Tower Surge Impedance," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-88, No. 5, pp. 680-687, 1969. [30] A . Semlyen and A . Dabuleanu, "Fast and Accurate Switching Transient Calculations on Transmission Lines with Ground Return using Recursive Convolutions," IEEE Transactions on Power Apparatus and Systems, Vol.PAS-94, pp. 561-571, March/April 1975. [31] H.W. Dommel, "Line Model in EMTP," EMTP Review, October 1988 & January 1989. [32], J.R. Marti and H.W. Dommel, "Line Models for Lightning Studies," E L E C T R I C A L ASSOCIATION, March 1989. CANADIAN [33] J.R. Marti, et al, "Approximate Transformation Matrices for Unbalanced Transmission Lines," 9th Power System Conference, Lisbon 1987. [34] J.R. Marti, "Accurate Modelling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulations," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-101, pp 147-157, 1982. [35] L . Marti, "Simulation of Transients in Underground Cables with Frequency-Dependent Modal Transformation Matrices," IEEE Trans, on Power Delivery, Vol. 3, No. 3, pp. 109-1110, 1988. [36] H . W. Dommel, EMTP T H E O R Y BOOK, 2nd edn., Microtran Power System Analysis Corporation, Vancouver, B.C., Section 4, 1992. . [37] L . M . Wedepohl, H. V . Nguyen and G. W. Irwin, "Frequency-Dependent Transformation Matrices for Untransposed Transmission Lines using Newton-Raphson Method," Paper 95SM602-3 PWRS, IEEE PES Summer Meeting, Portland 1995. [38] G. Angelidis, A . Semlyen, "Direct Phase-Domain Calculation of Transmission Line Transients Using Two-Sided Recursions," IEEE Transactions on Power Delivery, Vol. 10, No. 2, pp. 941-949, April 1995. [39] T. Noda, N . Nagaoka and A . Ametani, "Phase Domain Modelling of Frequency-Dependent Transmission Lines by Means of an A R M A Model," Paper 95WM245-1 PWRD, IEEE PES Winter Meeting, New York 1995. [40] F. Castellanos, J. R. Marti, "Phase-domain Multiphase Transmission Line Modelling," IPST'95 - International Conference on Power System Transients, Lisbon, September 3-7, 1995, pp. 17-22. 133 [41] L . M . Wedepohl, Frequency Domain Analysis of Wave Propagation in Multiconductor Transmission Systems, E E 552 Class Notes and Computer Program, Dept. of Elec. Engr., The University of British Columbia, 1993. [42] V . Larsen of EFI in Trondheim, Norway; personal communication with H.W. Dommel. [43] Personal communication with H.W. Dommel. [44] The Fast Front Transient Task Force of the IEEE Modelling and Analysis of System Transients Working Group, "Modelling Guidelines for Fast Front Transients," Drafted version undergoing review process, No. FFT-01/02/95-40. [45] L . Marti, "Simulation of Electromagnetic Transients in Underground Cables with Frequency-Dependent Modal Transformation Matrices," Ph.D. dissertation, The University of British Columbia, November 1986. [46] H.V. Nguyen, H.W. Dommel, J.R. Marti, "Modelling of Single-Phase Nonuniform Transmission Lines in Electromagnetic Transient Simulations," Accepted for presentation, 1996 IEEE Summer Power Meeting, Denver, Colorado. [47] H.V. Nguyen, H.W. Dommel, J.R. Marti, "Direct Phase-Domain Modelling of Frequency-Dependent Overhead Transmission Lines," Accepted for presentation, 1996 IEEE Summer Power Meeting, Denver, Colorado. [48] Personal communication with L . M . Wedepohl. [49] L . M . Wedepohl, "Transient Analysis of underground power-transmission systems System-model and wave-propagation characteristics," Proc. IEE, vol. 120, No. 2, pp. 253-260,1973. [50] J.R. Carson, "Wave propagation in overhead wires with ground return," Bell Syst. Tech. Journal, vol. 5, pp. 539-554, 1926. [51] A. Deri et al, "A Simplified Model for Homogeneous and Multi-Layer Earth Return," IEEE Trans, on Power Apparatus and Systems, Vol. PAS-100, No. 8, pp. 3686-3693, 1981. [52] C. Dubanton, "Calcul approche des parametres primaires et secondaraires d'une ligne de transport. Valeurs homopolaires" ("Simplified calculation of primary and secondary parameters of transmission line. Zero sequence values", in French), E.D.F. Bulletin de la Direction des Etudes et Recherches, Serie B, pp. 52-62, 1969. [53] L . M . Wedepohl, "Application of matrix methods to the solution of travelling-wave phenomena in polyphase systems," Proc. IEE, vol. 110, pp. 2200-2212, Dec. 1963. 134 [54] R.H. Galloway, W.B. Shorrocks and L . M . Wedepohl, "Calculation of electrical parameters for short and long polyphase transmission lines," Proc. IEE, vol. I l l , pp. 2051-2059, Dec. 1964. [55] H . W. Dommel, Fourier Transformation Program, Microtran Power System Analysis Corporation, Vancouver, B.C. [56] H.V. Nguyen, Discussion of "Computation of Electromagnetic Transients on Nonuniform Transmission Lines," Paper 95 S M 397-0-PWRD, IEEE PES Summer Meeting, Portland, 1995. [57] A.J.F. Keri, Y.I. Musa, and J.A. Halladay, "Insulation Coordination for Delta Connected Transformers," IEEE Trans, on Power Delivery, Vol. 9, No. 2, pp. 772-780, April 1994.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Simulation of lightning surges on transmission lines
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Simulation of lightning surges on transmission lines Nguyen, Huyen V. 1995
pdf
Page Metadata
Item Metadata
Title | Simulation of lightning surges on transmission lines |
Creator |
Nguyen, Huyen V. |
Date Issued | 1995 |
Description | Computer simulations with programs such as the EMTP are widely accepted for lightning transient studies, to design the system in such a way that equipment failures and transmission line outages are prevented. In these studies, an accurate representation of the transmission lines and towers is required. Two new and improved time-domain models for the representation of these components are presented in this thesis. They are the frequency- dependent and the exponential transmission line models. The former is based on synthesizing the line functions directly in phase co-ordinates, and is used to represent multiphase transmission lines, while the latter is a single-phase model, which also synthesizes the line functions, to be used for tower representations. These models are valid for the required frequency range of 1 Hz to 1 MHz. They can be used for lightning studies as well as for many other types of transient studies. Extensive time-domain simulations have been performed to compare the results with those produced by existing line models. The results indicate that the new models are more accurate than some of the existing models, and numerically stable. Simulations of typical lightning surges have also been performed to compare the solutions produced by the developed models with actual measurements. With these new models, simulations match field recorded transient waveforms very well. In addition, a simple, fast and accurate method based on the Newton-Raphson iteration technique is introduced for finding transmission line modal parameters. It overcomes the eigenvector and eigenvalue switching problems experienced by most other methods. With a proper constraint equation for scaling the eigenvectors, this proposed method produces modal parameters which can reasonably be fitted with rational functions of minimum-phase type. With this, the modal transmission line modelling can be improved. The routine is very stable and can be used to diagonalize the Y Z product for uncommon lines with strong asymmetry. It has also been used to obtain the phase-domain parameters for the new frequency- dependent line model. |
Extent | 4968449 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-02-19 |
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.0065121 |
URI | http://hdl.handle.net/2429/4816 |
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 |
Graduation Date | 1996-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1996-091414.pdf [ 4.74MB ]
- Metadata
- JSON: 831-1.0065121.json
- JSON-LD: 831-1.0065121-ld.json
- RDF/XML (Pretty): 831-1.0065121-rdf.xml
- RDF/JSON: 831-1.0065121-rdf.json
- Turtle: 831-1.0065121-turtle.txt
- N-Triples: 831-1.0065121-rdf-ntriples.txt
- Original Record: 831-1.0065121-source.json
- Full Text
- 831-1.0065121-fulltext.txt
- Citation
- 831-1.0065121.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065121/manifest