A Frequency-dependent Multiconductor TransmissionLine Model with Collocated Voltage and CurrentPropagationbyArash TavighiM.Sc. University of Tehran, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Electrical and Computer Engineering)The University of British Columbia(Vancouver)March 2017© Arash Tavighi, 2017AbstractThis research contributes to developing a time domain and a frequency domain for-mulations to solve electromagnetic transients in power system with multiconductoroverhead transmission lines.The time domain solution introduces a frequency dependent transmission linemodel “FDLM”. For the development of the FDLM a fundamental constraint isadded to the classical line equations to maintain the symmetry between electric andmagnetic fields. As a result, voltage waves and current waves travel together andthe characteristic impedance remains uniform along the line. With this premise, aconstant real transformation matrix can be obtained to diagonalize the line func-tions with high accuracy. This feature can greatly facilitate the line modelling asopposed to the existing line models which require complex frequency dependenttransformation matrices for their diagonalization. The use of a single constant realtransformation matrix for the voltage and current waves which is exact over the fre-quency range enables FDLM to provide higher accuracy and numerical efficiencythan the existing line models while it complies with the physical system.The accuracy of the FDLM is assessed through comparisons with a newly de-veloped Discrete Time Fourier Series frequency domain solution. This methodol-ogy is based on the correct specification of the time window and frequency windowwidths. Guidelines are provided for this set up which avoids the typical Gibbs andaliasing errors related to the classical frequency domain solutions. The proposedfrequency domain solution is simpler to implement than the most commonly usednumerical Laplace transform solution while it does not require further considera-tions to use damping factors or windowing functions.iiPrefaceBased on the research presented in this dissertation, several papers have been pub-lished in and/or submitted to conference proceedings and journal articles. In allof the publications, I have implemented the models, developed the case studies,conducted simulations, analyzed results, prepared the initial manuscript (exceptfor the last paper). My supervisor Prof. Jose´ R. Martı´ has provided me instructivecomments and corrections throughout the process of conducting research studies,preparing and editing manuscripts.Dr. J. A. Gutierrez-Robles from the University of Guadalajara provided theLaplace method for our benchmarking purposes. Vero´nica A. Galva´n, a visitingscholar from Cinvestav, extended the Laplace method for various case studies. Pro-fessor Hermann W. Dommel provided us many instructive discussion sessions ontransients studies and transmission line modelling.The following describes published and submitted papers.- Chapter 2 was presented at the IPST 2015 conference. A. Tavighi, J. R.Martı´, and J. A. Gutierrez-Robles, “Comparison of the fdLine and ULM FrequencyDependent EMTP Line Models with a Reference Laplace Solution”, InternationalConference on Power Systems Transients (IPST2015), Cavtat, Croatia, June 15–18,2015.- Part of Chapter 3 submitted as a journal article. A. Tavighi, V. A. Galva´n, J. R.Martı´, and J. A. Gutierrez-Robles “Validation and Usage Guidelines for DiscreteTime Fourier Series to Solve Electromagnetic Transients in Power Systems”.- Part of Chapter 3 submitted as a conference paper. A. Tavighi, J. R. Martı´, V.A. Galva´n, J. A. Gutierrez-Robles, and H. W. Dommel, “A Discrete Time FourierSeries Formulation to Simulate Electromagnetic Transients for Multi-conductoriiiTransmission Lines”.- Chapter 4 submitted as a journal article. J. R. Martı´, and A. Tavighi, “Fre-quency Dependent Multiconductor Transmission Line Model with Collocated Volt-age and Current Propagation”.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 History of Line Parameters Calculation . . . . . . . . . . . . . . 31.3 History of Frequency Domain Solutions . . . . . . . . . . . . . . 61.4 History of Time Domain Line Models . . . . . . . . . . . . . . . 71.5 Research Objectives and Anticipated Impacts . . . . . . . . . . . 111.5.1 Validation of Frequency-Dependent Line (FDLINE) UnderAsymmetrical Line Configurations . . . . . . . . . . . . . 121.5.2 Development of an Accurate Frequency Domain Solution 12v1.5.3 Revisions to the Classical Multiconductor Transmission Line(MTL) Equations and Development of a Frequency Depen-dent Electro-Magnetic Transients Program (EMTP) Line ModelBased on the Revised Equations . . . . . . . . . . . . . . 122 Accuracy of FDLINE Under Asymmetrical Line Configurations . . . 142.1 Classical MTL Equations . . . . . . . . . . . . . . . . . . . . . . 152.1.1 Per Unit Length Z and Y Calculation . . . . . . . . . . . 152.1.2 Diagonalization of the MTL Equations . . . . . . . . . . . 192.2 Time Domain and Frequency Domain Solutions to the MTL Equations 212.2.1 FDLINE Solution to the MTL Equations . . . . . . . . . . 212.2.2 Universal Line Model (ULM) Solution to the MTL Equations 242.2.3 Numerical Laplace Transform (NLT) Solution to the MTLEquations . . . . . . . . . . . . . . . . . . . . . . . . . . 252.3 Case Studies and Simulation Conditions . . . . . . . . . . . . . . 282.4 Comparison of Simulation Results . . . . . . . . . . . . . . . . . 332.4.1 Single-Phase Line . . . . . . . . . . . . . . . . . . . . . 332.4.2 Three-Phase Single-Circuit Lines . . . . . . . . . . . . . 342.4.3 Three-Phase Double-Circuit Line in the Same Tower . . . 352.4.4 Three-Phase Double-Circuit Line in Separate Towers . . . 372.4.5 Steady-state Solutions . . . . . . . . . . . . . . . . . . . 382.5 General Observations . . . . . . . . . . . . . . . . . . . . . . . . 383 Discrete Time Fourier Series Algorithm . . . . . . . . . . . . . . . . 403.1 Relationship Between Discrete Time Domain and Discrete Fre-quency Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.2 Power System Transient Solution Using the Discrete Time FourierSeries (DTFS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.2.1 Defining the Time Window Width (Tc) . . . . . . . . . . . 453.2.2 Defining the Frequency Window Width ( fc) . . . . . . . . 453.2.3 System Solution Using the DTFS . . . . . . . . . . . . . . 463.3 Guidelines to Specify Tc and fc . . . . . . . . . . . . . . . . . . . 473.3.1 Calculation of the Time Window Width (Tc) . . . . . . . . 47vi3.3.2 Calculation of the Frequency Window Width ( fc) . . . . . 483.4 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.4.1 Single-phase Circuits . . . . . . . . . . . . . . . . . . . . 503.4.2 MTL-Circuits . . . . . . . . . . . . . . . . . . . . . . . . 513.5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . 543.5.1 RL Circuit Simulation . . . . . . . . . . . . . . . . . . . 573.5.2 Nominal-pi Circuit Simulation . . . . . . . . . . . . . . . 623.5.3 Exact-pi Circuit Simulation . . . . . . . . . . . . . . . . . 623.5.4 MTL-Circuits Simulation . . . . . . . . . . . . . . . . . . 643.6 Sensitivity of the DTFS to the Width of the tset . . . . . . . . . . . 703.7 General Observations . . . . . . . . . . . . . . . . . . . . . . . . 723.8 Sensitivity of FDLINE to the Frequency of Approximation for theTransformation Matrices . . . . . . . . . . . . . . . . . . . . . . 724 Revised Multiconductor Transmission Line Equations, Proposed Frequency-Dependent Line Model . . . . . . . . . . . . . . . . . . . . . . . . . 784.1 Inconsistencies in the Classical MTL Equations . . . . . . . . . . 794.2 Collocation of Current and Voltage Waves . . . . . . . . . . . . . 804.3 Effect of the Earth Return in the Matrix of Impedances . . . . . . 814.4 Validity of the Collocation Condition ZY = YZ . . . . . . . . . . 844.5 Diagonalization of the Revised Multiconductor Transmission Line(RMTL) Equations . . . . . . . . . . . . . . . . . . . . . . . . . . 854.6 Transformation Matrix T to Diagonalize Z and Y . . . . . . . . . 864.7 Transformation Matrix T to Diagonalize ZY and YZ . . . . . . . 894.8 Proposed Frequency-Dependent Line Model . . . . . . . . . . . . 914.9 Numerical Results: Transient Simulations . . . . . . . . . . . . . 944.10 General Observations . . . . . . . . . . . . . . . . . . . . . . . . 984.11 Further Investigations on the Wave Collocation Condition . . . . . 1005 Summary of Contributions and Future Works . . . . . . . . . . . . 1075.1 Conclusions and Contributions . . . . . . . . . . . . . . . . . . . 1075.2 Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111viiList of TablesTable 2.1 Comparison of steady state values of different solutions. . . . . 39Table 3.1 Analytical formulas for the RL and nominal-pi circuits. . . . . 52Table 3.2 Modal lumped elements of the nominal-pi circuit for the doublecircuit vertical line of Fig. 3.5 calculated at 10−4 Hz. . . . . . . 52Table 3.3 Time-frequency parameters for the DTFS. . . . . . . . . . . . . 56Table 4.1 Errors in the diagonalization of the matrix of ones [1] and Cextfor T calculated at 100 Hz. . . . . . . . . . . . . . . . . . . . 89Table 4.2 Real part of the eigenvalues of ZY calculated using the pro-posed single T versus exact diagonalization (all values are neg-ative). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Table 4.3 Imaginary part of the eigenvalues of ZY calculated using theproposed single T versus exact diagonalization (all values arenegative). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Table 4.4 Maximum errors for voltages and currents for FDLINE, ULM,and Frequency-Dependent Line Model (FDLM) with respect totheir reference DTFS solution. . . . . . . . . . . . . . . . . . . 96Table 4.5 Comparison of operations count for FDLINE, ULM, and FDLMline models in the simulation of unbalanced faults test of Fig. 3.6. 97viiiList of FiguresFigure 2.1 Transmission line segment ∆x (only one conductor with re-spect to ground is shown). . . . . . . . . . . . . . . . . . . . 15Figure 2.2 Conductors above non-ideal ground using the method of im-ages and complex penetration depth. . . . . . . . . . . . . . . 17Figure 2.3 FDLINE modal equivalent circuit. . . . . . . . . . . . . . . . 22Figure 2.4 Coupled exact-pi model. Frequency domain equivalent-circuitfor multiconductor transmission line. . . . . . . . . . . . . . 27Figure 2.5 Single-phase transmission line (open and shorted). . . . . . . 29Figure 2.6 Three-phase single-circuit horizontal transmission line. . . . . 29Figure 2.7 Three-phase single-circuit vertical transmission line. . . . . . 30Figure 2.8 Three-phase single-circuit delta transmission line. . . . . . . . 30Figure 2.9 Three-phase double-circuit one-tower delta transmission line. 31Figure 2.10 Three-phase double-circuit two-tower horizontal transmissionline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Figure 2.11 One-line diagram of the equivalent circuit for the tests. a)Single-circuit, b) Double-circuit. . . . . . . . . . . . . . . . . 32Figure 2.12 Simulation results for the single-phase line of Fig. 2.5. . . . . 33Figure 2.13 Simulation results for the horizontal line of Fig. 2.6. . . . . . 34Figure 2.14 Simulation results for the vertical line of Fig. 2.7. . . . . . . . 35Figure 2.15 Simulation results for the delta line of Fig. 2.8. . . . . . . . . 36Figure 2.16 Simulation results for the double circuit delta line of Fig. 2.9. 37Figure 2.17 Simulation results for the two parallel horizontal line of Fig.2.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38ixFigure 3.1 Correspondence of time domain and frequency domain for adiscrete signal. a) Ns = even. b) Ns = odd. . . . . . . . . . . . 42Figure 3.2 Time window selection. a) System impulse response h[n], b)Current source i[n] (Dashed line: unit step with discontinuity,Solid line: unit step with the averaging of points at the discon-tinuity). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Figure 3.3 Example single-phase overhead transmission line. . . . . . . . 50Figure 3.4 Single-phase line models. a) exact-pi , b) nominal-pi , c) RL. . . 50Figure 3.5 Double circuit vertical transmission line. . . . . . . . . . . . . 53Figure 3.6 System diagram for the unbalanced faults test. . . . . . . . . . 54Figure 3.7 System diagram for the induced voltage test. . . . . . . . . . 54Figure 3.8 System diagram for the step response test. . . . . . . . . . . . 55Figure 3.9 isc for the RL circuit energized with a step voltage. . . . . . . 57Figure 3.10 Close up of isc for the RL circuit energized with the step volt-age for tsim = 0.5 s. . . . . . . . . . . . . . . . . . . . . . . . 58Figure 3.11 Effect of different windowing functions to reduce the oscilla-tions of the NLT solution in Fig. 3.10. . . . . . . . . . . . . . 59Figure 3.12 vL for the RL circuit energized with the step voltage. . . . . . 59Figure 3.13 Test of NLT with odd sampling in the simulation of vL for theRL circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Figure 3.14 The effect of filtering on the accuracy of NLT with odd sam-pling. Simulation of vL for the RL circuit. . . . . . . . . . . . 61Figure 3.15 isc for the nominal-pi circuit energized with the step voltage. . 62Figure 3.16 Details of Fig. 3.15 illustrating the effect of the first-last pointcorrection. . . . . . . . . . . . . . . . . . . . . . . . . . . . 63Figure 3.17 isc for the exact-pi circuit energized with the step voltage. . . . 63Figure 3.18 isc for the exact-pi circuit energized with the cosine voltage. . . 64Figure 3.19 Error in the DTFS due to incorrect selection of the time window. 65Figure 3.20 v3 in the unbalanced faults test of Fig. 3.6. . . . . . . . . . . . 66Figure 3.21 i6 in the unbalanced faults test of Fig. 3.6. . . . . . . . . . . . 67Figure 3.22 v6 and i2 in the induced voltage test of Fig. 3.7. . . . . . . . . 67Figure 3.23 v1 in the step response test of Fig. 3.8. . . . . . . . . . . . . . 68xFigure 3.24 The real part of TI(ω) for the double circuit vertical line ofFig. 3.5. a) with modal switching. b) corrected for modalswitching. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 3.25 v4 in the step response test of Fig. 3.8 illustrating the effect ofmodal switching problem. . . . . . . . . . . . . . . . . . . . 69Figure 3.26 Sensitivity of the DTFS method to the width of the tset for dif-ferent tests. a) Relative error, b) Computational time. . . . . . 71Figure 3.27 Sensitivity of FDLINE to the frequency of approximation . . . 74Figure 3.28 Close up of v3 and i6 in the unbalanced faults test of Fig. 3.6. . 75Figure 3.29 Close up of i2 and v6 in the induced voltage test of Fig. 3.7. . 76Figure 3.30 v1 in the step response test of Fig. 3.8. . . . . . . . . . . . . . 77Figure 4.1 Test circuit to determine the self and mutual impedances in theseries impedance matrix Z. . . . . . . . . . . . . . . . . . . . 82Figure 4.2 Comparison of Re in MTL and RMTL models for the doublecircuit vertical line of Fig. 3.5. . . . . . . . . . . . . . . . . . 84Figure 4.3 Validity of the collocation condition. If ZY = YZ, the productZY(YZ)−1 should be one for the diagonals and zero for theoff-diagonals. . . . . . . . . . . . . . . . . . . . . . . . . . . 85Figure 4.4 RMTL formulation. T matrix to diagonalize Lloop for the lineof Fig. 3.5. The same matrix T diagonalizes ZY and YZ. . . . 87Figure 4.5 MTL formulation. Real and imaginary part of complex TV andTI matrices to diagonalize ZY and YZ for the line of Fig. 3.5. 88Figure 4.6 Diagonalization of T−1ZYT and T−1YZT. The plots show theabsolute value of the off-diagonal elements. . . . . . . . . . . 90Figure 4.7 FDLM equivalent circuit in the frequency domain. Quantitiescan be matrices in the phase domain or scalar quantities in themodal domain. . . . . . . . . . . . . . . . . . . . . . . . . . 93Figure 4.8 Comparison of MTL (FDLINE, ULM) and RMTL (FDLM) linemodels for the unbalanced fault test of Fig. 3.6. . . . . . . . . 95Figure 4.9 Details of the first four peaks of v3 in Fig. 4.8a. . . . . . . . . 96xiFigure 4.10 Errors for the voltage v3 for FDLINE and ULM versus DTFS forthe MTL equations and for FDLM versus DTFS for the RMTLequations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Figure 4.11 v1 in the step response test of Fig. 3.8. Comparison of RMTLmodels: DTFS versus FDLM (for FDLM T calculated at differentfrequencies). . . . . . . . . . . . . . . . . . . . . . . . . . . 98Figure 4.12 Earth transversal losses given by the complex penetration depthformula for the double circuit vertical line of Fig. 3.5 . . . . . 101Figure 4.13 Lloop and Cloop for the line of Fig. 3.5. . . . . . . . . . . . . 102Figure 4.14 Validity of the collocation condition. If ZY = YZ, the productZY(YZ)−1 should be one for the diagonals and zero for theoff-diagonals. . . . . . . . . . . . . . . . . . . . . . . . . . . 103Figure 4.15 The effect of skin effect in capacitances on the product LC forthe double circuit vertical line of Fig. 3.5. . . . . . . . . . . . 104Figure 4.16 Comparison of MTL and RMTL models for two conditions: a)no skin in C, b) with skin in C. v3 in the test of Fig. 3.6. . . . 105Figure 4.17 Comparison of MTL with Cext versus RMTL with Cloop. v3 inthe test of Fig. 3.6. . . . . . . . . . . . . . . . . . . . . . . . 105xiiGlossaryACSR Aluminum Conductor Steel ReinforcedBAF Bode’s Asymptotic FittingBPA Bonneville Power AdministrationDTFS Discrete Time Fourier SeriesEMTP Electro-Magnetic Transients ProgramFDLINE Frequency-Dependent LineFDLM Frequency-Dependent Line ModelFEM Finite Element MethodFFT Fast Fourier TransformIDTFS Inverse Discrete Time Fourier TransformIFFT Inverse Fast Fourier TransformMFT Modified Fourier TransformMTL Multiconductor Transmission LineNLT Numerical Laplace TransformNSERC Natural Science and Engineering Research CouncilPPC Points Per CyclexiiiRMTL Revised Multiconductor Transmission LineTEM Transverse Electric MagneticUBC University of British ColumbiaULM Universal Line ModelVF Vector-FittingxivAcknowledgmentsForemost, I would like to express my sincere gratitude to my supervisor Prof. Jose´R. Martı´ for his continuous support, patience, great inspiration, immense knowl-edge and insights, and enthusiasm that he provided me throughout my studies. Iappreciate all his contribution of time and ideas to make my doctoral experienceproductive and stimulating. I am also grateful to him for the knowledge that I ac-quired from his “Network Analysis and Simulation” course, and for the wonderfulteaching experience that I gained while assisting him for the “Power System Anal-ysis 1 and 2” courses. The financial support for this research was made possible bythe Natural Science and Engineering Research Council (NSERC) of Canada undera grant led by Prof. Jose´ R. Martı´ as a principal investigator.I would like to express my appreciation to Dr. J. A. Gutierrez-Robles andVero´nica A. Galva´n for their valuable discussions on frequency domain solutionsduring their visit at the University of British Columbia (UBC).My sincere thanks also go to Prof. Hermann W. Dommel for his valuablesuggestions and knowledge contributions during the development of this research,and to Prof. Juri Jatskevich for his constructive comments on the presentation ofcontents of this dissertation and his magnificent lectures on “Dynamic Modellingof Electrical Machines and Controls” course.Finally, I would like to thank my colleagues graduate students in Electric Powerand Energy Systems research group at UBC: William Wang, Pouya Zadkhast, Mo-hammad Ghasemi, Aurash Alimardani, Francis Therrien, Ehssan Ghahremani, andAndrea Martı´ for many stimulating discussions and all the memorable momentswe have had in the past five years at UBC.xvDedicationTo Maria and My Beloved ParentsxviChapter 1Introduction1.1 MotivationOverhead transmission lines and underground cables transfer electric power overmany kilometers from the generation sector to the customers at the distributionlevels. These links are usually affected by over-voltages and over-currents causedby wide variety of transients such as switching operations, short circuits, lightningdischarges, and faulty insulation. If these stresses imposed to the transmission linesare beyond a sustainable limit, these links can be severely damaged. To restrain theamplitude and duration of these undesirable phenomena, protective devices suchas surge arresters and fault detection relays should be employed, and their locationand setting should be tactfully determined. The success of the protective strategiesto disconnect the lines during the transients depends on the accurate prediction ofthe system. This can be done through simulations when accurate transmission linemodels are available.Two main approaches can be used to solve electric transients in frequency-dependent transmission lines: Time domain and frequency domain solutions.The frequency domain solution describes the transmission line one frequencyat the time over a frequency range, and transforms the frequency response into atime response using time-frequency transformations such as the Fourier Transform.These models can analytically describe the frequency dependence of the circuitparameters, and unlike the time domain solutions they make no approximations1when bridging the time response to the frequency response. Therefore, frequencydomain solutions are valuable tools to assess the accuracy of the correspondingtime-domain models. However, frequency domain solutions have difficulty to rep-resent directly the opening and closing of circuit breakers at specific times andmodel nonlinear elements. Simulation of switching operations in frequency do-main solutions requires going back and forth between frequency and time domainsfor each change of switch position, thus, making the process quite long and dif-ficult. In terms of the computational costs, frequency domain solutions are moreexpensive than the time domain solutions as they require a one-to-one correspon-dence between frequency and time domains. This restricts their applications to beused as general purpose solvers.Time domain methods are more flexible for circuits with switching operations,non-linear and time-varying elements, and to perform real-time analysis. Amongthe time-domain methods, the Electro-Magnetic Transients Program (EMTP) [29][30] [28] is the most widely used tool. In direct time domain simulations, everyelement of an electrical network is described by a set of differential equations.In order to obtain the time-domain models of these elements, the correspondingdifferential equations are solved (integrated) between discrete time intervals (e.g.trapezoidal rule in the EMTP. The main handicap of time-domain methods is mod-elling frequency dependent components. Nonetheless, very good frequency de-pendent transmission line and underground cable models exist in the EMTP, as forexample, Frequency-Dependent Line (FDLINE) [69], fqLine [71], and UniversalLine Model (ULM) [74]. These models approximate the frequency dependence ofthe line or cable wave functions with rational function approximations (using, forexample, Bode’s Asymptotic Fitting (BAF) [1] and Vector-Fitting (VF) [55]). Thefitted functions can then be transferred into the time domain to incorporate thesemodels into the time domain solution. This conversion has to be done numerically,as a result, selection of the time step and numerical stability can affect the accuracyof the modelling [68]. Besides other factors such as consideration of the transfor-mation matrices to diagonalize the line functions, the accuracy of the curve-fittingalgorithms can also determine the accuracy of these line models.The problem of frequency dependent transmission line modelling has been thesubject of many researches for many years. Despite the significant progress made2in this area, there is still room to further advance the modelling by taking intoaccount physical assumptions.In this section, the advances in the calculation of line parameters used in theclassical Multiconductor Transmission Line (MTL) equations are briefly reviewed.The history of the most well-known frequency domain solutions and time domainline models are presented along with their pros and cons in preparation of propos-ing the Revised Multiconductor Transmission Line (RMTL) equations as the mainobjective of this research.1.2 History of Line Parameters CalculationTraditional ApproachThe MTL equations belong to the class of formulations where the skin effect isonly considered in the series impedance; while the shunt admittance is calculatedassuming electrostatic position for the charges.Sommerfeld in 1909 proposed a pioneering solution to the problem of “elec-tromagnetic wave propagation along current-carrying conductors above a finitelyconducting plane” [112]. In 1926, Carson proposed Quasi-Transverse ElectricMagnetic (TEM) fields approximations of the Maxwell equations to solve this prob-lem for an overhead transmission line [13]; whereas Pollaczek developed similarformulas applicable not only to overhead transmission lines but also to the un-derground cables and to combinations of both [98] [97]. Although Carsons ap-proximations are relatively simple, they are valid for limited range of frequenciesonly that the capacitive displacement currents in the ground can be neglected. Tovalidate Carsons equations for high frequencies, Semlyen presented an analyticalcontinuation to the Carson’s integral [106], and Hofmann derived an asymptoticseries expansion for the calculation of the self and mutual line impedances [58].There have been several attempts to approximate Carsons integral aiming tolessen the computational time of the simulation. First, Dubanton in 1969 sug-gested a simple and sufficiently accurate expressions for line impedances for thewhole range of frequencies with intuitive insight [31]. Dubanton’s formula wasanalytically proven by Gary in 1976 [41] and Deri in 1981 [26] introduced as the“complex ground return plane approach” with further improvement in [105] [110]3[4]. Other authors presented approximations to Carson’s formula in a form ofdouble logarithmic [96] with further enhancement to obtain a simpler and moreaccurate closed-form formula in [87], and a more recent formulation based on thevariable changes integral in [61]. Numerical evaluation of the Carson’s integralwere presented in [100] [86] [64] to assess the accuracy of approximation formu-las.Stratified EarthCarson’s formula was developed for the homogeneous infinite earth. In prac-tice, however, the earth is not homogeneous, and its resistivity varies along thedepth of the earth layer. Even if considered as an equivalent homogeneous earth,the resultant earth resistivity may be frequency dependent.In 1949, Sunde [113] formulated the series impedance to include a two-layerearth only for overhead lines or cables above the earth. In his solution, however,the boundary conditions are not sufficiently general, and the propagation of currentalong a line is neglected [83]. The idea of Carson’s integral extension for a two-layer earth was followed in [59] [78]. In 1973, Nakagawa proposed a more rigoroussolution for a three layers earth case [83], with a sensitivity analysis for arbitraryearth resistivities, pemeabilities and permittivities [6], and with further extensionto a general multi-layer case with continuously variable (exponential or linear)earth resistivity [82]. Nakagawa found out that the homogeneity assumption isonly permissible at very high frequencies. Indeed, at low frequencies, a stratifiedearth causes significant differences in the earth impedances and the resultant wavedeformations from the homogeneous case. He also concluded that the displacementcurrents mostly affect the earth-return impedances at frequencies beyond 1 MHzand under the conditions that the earth resistivity is high and the height of theconductors above the ground is low [83].Further studies in the modelling of the stratified earth include a direct numericalintegration using Finite Element Method (FEM) [91] and closed-form formulas fora two layers earth [8] [119] with further generalization for a multi-layer earth [121][120].Correction on AdmittanceAs discussed, the MTL equations are based on the TEM assumption under whichthe ground return correction is only applied to the series impedance; while the shunt4admittance is considered only of the physical geometry. Under this condition, it isassumed that the conductor and the earth surface have the same potential [40].As a result, only at high frequencies where the correction term of the impedancebecomes negligible, the admittance agrees with the impedance and the propagationconstant becomes purely imaginary and waves propagates with no loss above lossyground [80]. However, when the ground is not perfectly conductive (non-idealground), the potential at the surface of the ground should not be considered aszero. Therefore, ground return corrections have to also be considered for shuntadmittances of transmission lines [81].In 1930’s, Wise extended Carson’s approach for relative permeability and per-mittivity of the earth greater than the unity [139] [137], and added the effects oftransversal currents displacement for a single-conductor above an imperfect earth[138]. He referred to restricting assumptions of Carson’s approach such as neglect-ing the polarization currents at low frequencies, and propagation of the wave withthe velocity of the light and without attenuation. Kikuchi is the first author whoproved the exact theory to Wise’s solution by calculating the admittance correctionterm by means of a series expansion and pointing out the surface wave propagationat high frequencies [62]. Kikuchi’s formula agrees with Wise’s under the conditionthat the permeabilities of both medias are equal.Sunde [113] and Arisumunandar [9] proposed similar image approximation[26] for the longitudinal effects of a transmission line. Hedman performed modalanalysis on the admittance correction terms of Wise for high frequencies [56]. Heconcluded that the effects of the high relative dielectric constant of the earth aresignificant only for frequencies higher than 0.5 MHz, and when both earth resis-tivity and dielectric constant are high. Wedepohl further extended Wise correctionterm for a semi-infinite non-homogeneous earth [133] in terms of an infinite inte-gral with no restriction on the permeabilities, permittivities, and resistivities. Hisnumerical results showed that the integral becomes equivalent to that derived byCarson and Wise for the case of a semi-infinite homogeneous earth. However, con-siderable differences have been shown to occur in the stratified cases, when theearth-layer resistivities differ by a factor of 10. Moreover, the capacitive nature ofthe earth was shown to be significant at very high frequencies due to displacementcurrents.5In 1969, Wait proposed an exact solution of an implicit modal equation basedon the full-wave approach of Sommerfeld integrals and Bessel functions in orderto leave out of the Carson approximation for a thin, infinitely long wire parallelto an interface [125]. The formulation is rather general but very difficult to han-dle from a computational point of view. D’Amore presented an approximation tothe exact solution of Wait for a single [21] and multi-conductor [22] transmissionlines with further development to obtain simpler formulation [20]. Wedepohl andEfthymiadis developed a more general numerical technique to evaluate propaga-tion constants for a single, lossless, infinitely long conductor above a lossy ground[129] [33]. This approach combines the vector potential as the sum of comple-mentary transverse magnetic and transverse electric wave components for low andhigh frequencies showing similar results to Kikuchi’s technique. The most obviousnew feature of the results is the increase in the propagation loss due to a changein the reactive nature of the ground when high ground resistivities were involved.They realized that introducing a negative conductance term which is normally neg-ligible, leads to a lossless propagation wave. This work was further enhanced byNakagawa for a single [81] and multi-conductor [80] cases resulting in an approachwith much smaller computation time than the solution of Efthymiadis. Nakagawarealized that in the case of ordinary lossy ground, admittance correction term doesnot have a significant impact in analyzing wave propagation at frequencies of in-terest to power engineers; whereas it produces remarkable changes in attenuationconstants at high frequencies and for very lossy medium.Pettersson generalized Sunde’s image representation approach for the trans-verse effects, and achieved simple closed-form expressions for the axial propaga-tion constant and characteristic impedance [94]. He further extended his formula-tion to include the vertical electric field between ground and wire [95]. Recently,Papadopoulos generalized Sunde’s admittance correction approximation for over-head conductors above a two-layer earth [90].1.3 History of Frequency Domain SolutionsAs mentioned earlier, frequency domain solutions are used to assess the accuracyof the time domain line models since they do not approximate the frequency de-6pendence of the line functions.Pioneering work on power system transients simulation using numerical Fouriertransforms was performed in 1965 [23]. This work addressed the Gibbs oscillationsdue to truncation of the integration range with infinite frequency of a continuous-time function. The authors introduced the use of the Lanczos’ window to reduceGibbs errors. They also referred to the aliasing errors in the time domain due tothe discretization of the continuous variable in frequency. To reduce time aliasing,they proposed the Modified Fourier Transform (MFT) in which they smoothenedthe frequency response by adding a damping factor [24]. Wedepohl applied modalanalysis to the MFT in the calculation of transients on multiconductor transmis-sion lines [132] and later for underground and submarine cables [134]. The resultsobtained by Wedepohl at that time could not be attained with the available time do-main methods and were used as benchmarks to enhance time domain line and cablemodels [101]. The computational efficiency of the MFT was enhanced by introduc-ing the application of Fast Fourier Transform (FFT) [5] and by using non-uniformsampling [7] [60] [18] [45].In 1978, Wilcox formulated the MFT methods in terms of the Laplace trans-form theory and introduced the term Numerical Laplace Transform (NLT) [136].He also provided a criterion to calculate the associated damping factor. In 1982,Wedepohl further enhanced the numerical accuracy of the NLT by calculating thedamping factor based on the number of samples [127]. The use of different win-dow functions and the choice of damping factors for the NLT were studied in [75].This reference also verified that the damping factor proposed by Wedepohl givesmore accurate results than the one proposed by Wilcox, at the cost of amplifyingthe Gibbs oscillations and numerical errors in the tails of the waveforms.Advances to include switching maneuvers, nonlinear and time-varying ele-ments in the frequency domain solutions can be found in [131] [130] [79] [77][122] [104] [76] [44] [109] [43] [42] [89] [17].1.4 History of Time Domain Line ModelsTime domain transmission line models are also referred to as “traveling-wave mod-els” in which the voltages and currents at the two line ends are decomposed into7incident and reflected waves. The behaviour of the line is described by the char-acteristic impedance which relates current waves to voltage waves and the propa-gation function which defines the delay and distortion of a wave traveling betweenthe two line ends. Traveling-wave line models can be also classified into threetypes: Constant Parameters Line (CP-Line) model, modal domain models, andphase domain models. The CP-Line model was initially introduced by Bergeronin 1961 [10], and was the first model incorporated in the EMTP. The CP-Line is asimple, single-line constant-frequency model. Here, the line is treated as losslesswith a pure time delay and a real characteristic impedance in the form of a dualThevenin equivalent circuits. However, its distributed series resistance is added inlump form. The lumped resistances can be inserted throughout the line by dividingits total length into several sections. The main disadvantage of the CP-Line is thatthe relatively strong frequency dependence of the line parameters is not taken intoaccount. Therefore, the CP-Line cannot adequately simulate the response of theline over the wide range of frequencies that are present in the signals during mosttransient conditions.The theory of modal decomposition was introduced by Wedepohl in 1963 [128]and Hedman in 1965 [56] (independently from each other). This theory is basedon the use of eigenvectors transformation matrices to decouple the physical systemof an N-phase multiconductor transmission line into N-mathematically-equivalentdecoupled single-phase circuits. Based on the modal decomposition theory and theCP-Line model, one of the first frequency dependent multiconductor line modelswas proposed by Budner in 1970 [12]. He introduced the concept of weightingfunctions to model the frequency dependence of the line characteristic admittanceand wave propagation functions. The weighting functions in this model are; how-ever, highly oscillatory and difficult to evaluate with accuracy [68]. In an effort toimprove Budner’s model, in 1972 Snelson [111] introduced an analogous variablechange to relate currents and voltages in the time domain. This idea was furtherdeveloped by Meyer and Dommel [73] in 1974 and resulted in obtaining forwardand backward travelling function from the weighted past history of the currentsand voltages at both line ends using the convolution integral. This formulation rep-resented a considerable improvement over other weighting function methods, andgave reliable results in many cases of transient studies [68]. On the other hand, the8relatively time consuming process required to evaluate the convolution integralsat each time step of the solution and the less accurate response at low frequen-cies, including the nominal 60 Hz steady state, were the main drawbacks of thistechnique [68]. As a solution, in 1975 Semlyen and Dabuleanu [107] introducedan efficient formulation to synthesize the line functions with a low-order rationalfunction approximation using complex exponentials. However, by allowing com-plex poles and zeros in the fitting functions, numerical oscillations were introducedin the response of line model.J. Martı´ in 1981 [68] [69] introduced the FDLINE model that simplified fre-quency dependence modelling by fitting the amplitude of line functions with highorder BAF rational approximations using negative poles and zeroes. This choiceresulted in a minimum-phase realization (i.e., phase function is uniquely deter-mined from the magnitude function) and, therefore, guarantees causality (i.e., thetransfer function is equal to zero for t < 0) and absolute numerical stability. InFDLINE, eigenvector transformation matrices are used to transform the originalcoupled phase coordinates into decoupled modal coordinates. FDLINE approxi-mates the complex frequency dependent transformation matrices obtained with theclassical MTL equations at one frequency (e.g., 1 kHz), and only takes the real partof the matrices. The original version of FDLINE was implemented in 1981-1983in the Bonneville Power Administration (BPA) EMTP. This version was improvedin the DCG/EPRI version developed in 1984-1986 [1]. Some of the specific im-provements included: higher dynamics in the very low frequency region that solvedencountered problems for very short lines, improved BAF algorithms that enhancedthe overall accuracy of the fitting, and automatic error checking for short/open cir-cuit conditions to optimize the single-frequency diagonalizing matrices used in thismodel. FDLINE is a simple, efficient, and reliable model that has been incorporatedin most versions of the EMTP (e.g., ATP, Microtran, PSCAD, EMTP-RV, and Opal-RT) and in real-time versions of the EMTP [32]. Even though FDLINE has beenwidely used, there has been a belief that because it uses a single real transforma-tion matrix to convert between modal and phase quantities, the model may not beaccurate enough for strongly asymmetrical line configurations [36] [38]. However,reference [70] showed that for switching surge studies in untransposed overheadlines, real constant transformation matrices can still give acceptable results within9certain frequency limits. A number of line models have been proposed to removethe limitation of FDLINE of using a real constant transformation matrix as they willbe briefly introduced next.L. Martı´ in 1988 [71] proposed the fqLine model which synthesized the trans-formation matrices with BAF rational function approximations with real negativepoles. This approach worked well in the case of underground cable where the con-ductors are very close to each other and to ground and the transformation matri-ces become strongly frequency dependent [72]. Fitting the transformation matriceswith minimum-phase rational functions proved a difficult task. Here, the numericaldifficulties of relaxing the constant transformation matrix condition for the over-head transmission lines began to emerge. Some of the terms that had to be fittedwith rational functions could not be satisfactorily approximated using only nega-tive poles and there was also the problem of modes switching along the frequencyrange [37].In 1999, Tavares used a transformation matrix based on line geometry andClarke transformation matrix to diagonalize frequency dependent multiphase trans-mission lines in the modal domain [114]. This model provides exact solution forideally transposed lines and single-circuit three-phase horizontal lines and providesgood approximation for non-transposed lines with vertical symmetry plane.Recently, Gustavsen showed that the accuracy of FDLINE in two parallel over-head lines can be improved by modeling each line by an FDLINE and then addingthe mutual coupling between the independent FDLINEs using rational functions[49]. However, this reference concluded that the proposed model seems to be com-putationally expensive particularly if the upper bandwidth passes 100 kHz.The difficulties involved in the modal domain line models led to the develop-ment of alternative line models that do not use transformation matrices by workingdirectly in phase coordinates. Examples of these models are the zLine [15], zCable[141], ARMA [88], [85], and other models based on the method of characteristics[84][102] [35].Other phase domain models are based on the idempotent theory which was in-troduced by Cullen in 1990 [19] as “spectral decomposition”. Wedepohl in 1993[126] used this concept to synthesize the propagation function of transmission linesin phase coordinates. The idempotents are slightly better behaved functions than10the eigenvectors of the modal transformation matrices; however, they are still fre-quency dependent. They also have the unique property to be scale-free, which is akey factor in solving the problem of normalization of the transformation matricesin fqLine model. Idempotent-based line models for the EMTP were first introducedby Castellanos and J. Martı´ in 1995 [14]. They were further developed by Marcanoand J. Martı´ [65] [66]. In these models accuracy was sacrificed in order to useonly minimum-phase functions for asymmetrical line configurations. Gustavsenand Semlyen in 1998 relaxed the requirement of minimum-phase in the synthesisof the idempotent coefficients to gain accuracy [52] [51] [54]. This work evolvedin the ULM [74] in 1999 which is widely used in the EMTP community when fre-quency dependence of the transformation matrices is a concern. In developing theULM, Gustavsen and Semlyen in 1999 introduced the VF rational approximationapproach [55] with the improvement in [46] [50]. By allowing unstable poles andzeros in the VF, unlike BAF in FDLINE, ULM requires both the real and imaginaryparts of the line functions to be fitted. Because of the higher burden of simulat-ing the frequency dependence of the transformation matrices, ULM is considerablymore expensive than FDLINE which is an important consideration in real time sim-ulators. VF requires post-processing to achieve passivity [53] [16] [48] [27] [25][108] and prevent numerical instability problems [63] [47]. Since physical systemswith no sources have to be passive, a mathematical model without intrinsic passiv-ity is not a good representation of reality and can lead to unexpected problems.1.5 Research Objectives and Anticipated ImpactsThe main objectives of this research are the followings:1. Validation of FDLINE under asymmetrical line configurations2. Development of an accurate frequency domain solution3. Revisions to the classical MTL equations and development of a frequencydependent EMTP line model based on the revised equations111.5.1 Validation of FDLINE Under Asymmetrical Line ConfigurationsThe FDLINE model has given reliable results in the EMTP community for manyyears. Even though FDLINE is simple, does not require passivity considerations,and has been implemented in real time simulators, the concern has remained ofwhether the model is accurate enough under asymmetrical line configurations.Many of the arguments for this conclusion have been based on using a single trans-formation matrix versus fitting the transformation matrix functions.In Chapter 2, the accuracy of the FDLINE is compared with the ULM and NLTfrequency domain solution for a variety of asymmetrical line configurations. Byvalidating the accuracy of FDLINE in this chapter, this model will be used in Chap-ter 4 in the development of a line model based on the RMTL equations.1.5.2 Development of an Accurate Frequency Domain SolutionIn Chapter 3, derivation of the Discrete Time Fourier Series (DTFS) methodologyis presented to solve electromagnetic transients in power systems with multicon-ductor lines. The proposed methodology is based on the correct specification ofthe time window and frequency window. Guidelines are provided for this set up.DTFS is a simpler solution than the most commonly used NLT due to not requiringto choose a damping factor and a windowing function to reduce aliasing and Gibbserrors typical of the frequency domain solutions; while it provides very accurateresults.1.5.3 Revisions to the Classical MTL Equations and Development of aFrequency Dependent EMTP Line Model Based on the RevisedEquationsChapter 4 addresses a physical inconsistency in the classical MTL equations thatis the voltage waves and current waves do not travel together along the line. As aresult, the characteristic impedance becomes a function of distance which contra-dictions Bergeron’s theory.To maintain the symmetry of propagation between the voltage waves and thecurrent waves, a physical constraint is imposed on the MTL equations. Based onthis condition the RMTL equations are proposed. As opposed to the MTL equationsthat require complex frequency dependent transformation matrices for their diag-12onalization, the RMTL equations can be diagonalized using a single real constanttransformation matrix.A Frequency-Dependent Line Model (FDLM) is proposed based on the RMTLequations. The accuracy of FDLM is compared to FDLINE, ULM, and DTFS fre-quency domain solution.13Chapter 2Accuracy of FDLINE UnderAsymmetrical LineConfigurationsThis chapter first reviews the classical MTL equations for an N-conductor over-head transmission line with ground return on which the existing line models arebased. This includes calculation of per unit length impedance and admittance ma-trices as well as derivation of the transformation matrices to diagonalize the wavefunctions. Next, well-known time domain and frequency domain solutions to theMTL equations are briefly discussed. Finally, the simulation results obtained withFDLINE and ULM line models in the EMTP are presented for a variety of asymmet-rical line configurations, including single-circuit, double-circuit and parallel lines.The results are compared with a conventional NLT reference solution. The com-parisons include open and short circuit conditions during the transient and steadystate periods.The MTL equations will be used in Chapter 3 in the development of the DTFSalgorithm, and in Chapter 4 to propose the RMTL equations and develop the FDLMmodel.142.1 Classical MTL Equations2.1.1 Per Unit Length Z and Y CalculationIn Fig. 2.1, a segment ∆x of a long uniform transmission line in a uniform mediumis shown. In general, the equations are derived for N parallel conductors above acommon ground return. The equations relating voltages and currents in this linecan be derived in the frequency domain in terms of series impedance and shuntadmittance matrices.IyxRcVyx∆LloopGins Cextx yZIx IyVx Ygx gyReFigure 2.1: Transmission line segment ∆x (only one conductor with respectto ground is shown).The series impedance and shunt admittance matrices in Fig. 2.1 are symmet-rical and their elements are complex numbers. Preserving this symmetry is oneof the physical constraints imposed in this research. Voltages are measured fromconductors to a single common ground. Currents in the conductors flow forwardin the x direction and return through the common ground path. The voltage andcurrent drops along the line conductors are given by,Vx−Vy = ZI , Ix− Iy = YV (2.1)The real part of the impedances in Z consists of the internal resistance of theconductors Rc plus the resistance of the ground return path Re. The imaginary partof Z consists of the inductances of the loops Lloop formed by the conductors andground. In the normal EMTP formulation, the real part of the admittances matrixY consists of the losses along the insulator strings Gins, and the imaginary part15consists of the external conductor capacitances Cext . The series impedance andshunt admittance matrices are then given byZ = (Rc+Re)+ jωLloop (2.2)Y = Gins+ jωCextIn general, all elements of the matrices in (2.2) can be frequency dependent.In normal EMTP transmission line modelling, however, Gins and Cext are assumedconstant.In (2.2), Rc and Gins are diagonal matrices. Rc is calculated using Bessel func-tions to account for the skin effect. In the EMTP simulation tools (e.g., Microtran,PSCAD, and EMTP-RV), Rc is calculated knowing the outer diameter, thicknessratio, and DC resistance of an Aluminum Conductor Steel Reinforced (ACSR) con-ductor. The value of the Gins is in the order of the (10−9)’s S/km [39] which isnormally small compared to the other system parameters. However, ignoring thisparameter causes a relative asymmetrical property between the Z and Y matrices,even though negligible.Re, Lloop, and Cext are full matrices which are calculated using the concept ofcomplex penetration depth [26].The diagram in Fig. 2.2 illustrates Lord Kelvin’s method of images wherefor ideal ground each line conductor at height hi has a corresponding image inthe ground at depth hi, and the actual ground can be replaced by the conductors’images.In an approximate solution, that illustrates the nature of the problem, the com-plex penetration depth method of [26] can be used to adjust the height of the con-ductors and the depth of their images (Fig. 2.2) to take into account non-idealground. The adjusting complex penetration depth value p is given byp=√ρjωµ(2.3)In (2.3), ρ is the earth resistivity, ω is 2pi f where f is the frequency and µ isthe permeability of the ground. For the air and the ground µ = µ0 = 4pi×10−7 inH/m. The radii of conductors is taken as rieq to correct for the internal skin effect.161Ground2d12h1 + p1'2'h1 + pD12'r1eqFigure 2.2: Conductors above non-ideal ground using the method of imagesand complex penetration depth.In Fig. 2.2 we show two conductors and their ground return images. In gen-eral we will have N physical conductors with N images. The classical MTL lineequations, however, are derived for N conductors above a single ground return.Mapping the system of conductors and images in Fig. 2.2 to the system of con-ductors with a common ground return, there is a loss of dimensionality of N− 1.This loss of information leads to the inconsistencies in the calculated values for Rewhich will be further discussed in Section 4.3.To calculate the values of Re using the method of images, we can define a com-plex loop inductance L′loop that takes into account the inductance and the resistanceof the path formed by the conductor and the earth return. The self element L′loop11in Fig. 2.2 is calculated inL′loop11 =µ2piln2(h1+ p)r1eq(2.4)The real part of the L′loop11 defines the self inductance of the loop, while theimaginary part corresponds to the self earth resistance given by17Lloop11 =ℜ(L′loop11)(2.5)Re11 = jω(jℑ(L′loop11))=−ωℑ(L′loop11) (2.6)The mutual element L′loop12 in Fig 2.2 is calculated inL′loop12 =µ2pilnD12′d12(2.7)In (2.7), D12′ is the distance from conductor 1 to the image of conductor 2, d12is the physical distance from conductor 1 to conductor 2.The mutual elements for inductance loop and earth resistance are given byLloop12 =ℜ(L′loop12)(2.8)Re12 = jω(jℑ(L′loop12))=−ωℑ(L′loop12) (2.9)For the Y matrix of (2.2), the matrix of capacitances Cext is obtained from thematrix of potential coefficients Pext asCext = P−1ext (2.10)Following the TEM condition in the classical MTL equations, the capacitanceis assumed to be equal to its electrostatic value (with the charges on the outsideof the conductor and on the surface of the earth). Since no transversal losses areconsidered, the values of Pext and Cext are real. The diagram of Fig. 2.2, with p=0, can be used to represent this electrostatic condition and calculate the potentialcoefficients. The self and mutual elements of real matrix Pext in the example of Fig2.2 are calculated inPext11 =12piεln2h1r1ext, Pext12 =12piεlnD12d12(2.11)In (2.11), r1ext is the external radius of conductor 1 in Fig. 2.2 and D12 is thedistance from conductor 1 to the image of conductor 2 when the ground penetrationdepth is zero. For the air and the ground ε = ε0 = 8.85×10−12 in F/m.182.1.2 Diagonalization of the MTL EquationsAnalogous formulas can be written for the voltage drops in the series impedancesbetween sections x and y and current drops across the shunt admittances, for con-ductors 1 and 2 in Fig. 2.1. Assuming Z is in Ω/m and Y is in S/m, and making thelength of the line segment ∆x go to zero, we obtain first order differential equationsfor the voltage and current drops as a function of the distance x down the line,− ∂V∂x= ZI , −∂ I∂x= YV (2.12)The negative signs in (2.12) are due to the assumed positive direction for thex in Fig. 2.1. Combining voltages and currents in (2.12), we obtain second orderdifferential equations for the propagation of voltage and current as a function ofthe distance x down the line,∂ 2V∂x2= (ZY)V ,∂ 2I∂x2= (YZ)I (2.13)Due to the duality of electric and magnetic field in the propagation process,symmetry conditions need to be preserved in the solution of the propagation equa-tions. Since the quantities in (2.13) are matrices, we need, in general, to preservethe order of the multiplications. It is then assumed in the classical MTL solution(e.g., [126][93]) thatZY 6= YZ (2.14)As we will show in Chapter 4, this assumption will lead to a physical con-tradiction in the symmetry of the solution. The solution of (2.13) can be writtenasV(x) = V f (x)+Vb(x) = V f ke−γVx+Vbke+γVx (2.15)I(x) = I f (x)+ Ib(x) = I f ke−γIx+ Ibke+γIx (2.16)The wave propagation constants γV and γI are found from the solutions of(2.13),γV =√ZY =√(Rc+Re+ jωLloop)(Gins+ jωCext) (2.17)19γI =√YZ =√(Gins+ jωCext)(Rc+Re+ jωLloop)(2.18)Without any further restrictions in the MTL equations, and except for the case ofa balanced line or a single-phase line, γV and γI in (2.17) and (2.18) are in generaldifferent.The voltage and current waves in (2.15) and (2.16) have two components: theforward waves “ f ” that travel from left to right in the diagram of Fig. 2.1, and thebackward waves “b” that travel from right to left in this diagram.Assuming a semi-infinite line, where only forward voltage and current wavesexist in (2.15) and (2.16), the voltage waves and the current waves are related by thecharacteristic impedance (2.19)-left. Similarly, the backward voltage and currentwaves are related as shown in (2.19)-right (the minus sign here indicates that thebackward currents flow from right to left),V f (x) = ZcI f (x) , Vb(x) =−ZcIb(x) (2.19)The characteristic impedance matrix Zc can be found by replacing (2.19) in(2.12) givingZc = Y−1√YZ (2.20)where Z and Y are defined in (2.2).For single-phase lines, Z and Y are scalar numbers and their product is commu-tative. In this case, γV and γI are the same. This is also the case for fully transposed(“balanced”) multiconductor lines that can be diagonalized by symmetrical com-ponents and other akin transformations. However, for untransposed multi-phaselines, γV and γI are different. This is the currently held assumption in the MTLequations.When modal analysis is used to transfer the line equations from phase to modalcoordinates (Wedepohl, 1963 [128] and Hedman 1965 [56]), two transformationmatrices TV and TI are defined, one for the voltages and one for the currents.Using subscript ph for the original coupled phase quantities and subscript m forthe decoupled modal quantities,20Vph = TVVm , Iph = TIIm (2.21)Matrix TV diagonalizes the product ZY while matrix TI diagonalizes the re-verse product YZ,ZYm = T−1V(ZphYph)TV , YZm = T−1I(YphZph)TI (2.22)In the modal domain, the voltage propagation ZYm is equal to the current prop-agation YZm.Since the transformation matrices contain complex numbers, they are prone tomode switching that has to be corrected. A well-known solution to this problemwas proposed by Wedepohl with the application of Newton-Raphson algorithm[135].The decoupled modal series impedances and shunt admittances are calculatedfromZm = T−1V ZphTI , Ym = T−1I YphTV (2.23)It will be argued in Chapter 4 that the constraint of collocation of voltage wavesand current waves will require γV and γI to be equal in the phase domain, regard-less of the asymmetry of the system of conductors (vertical, horizontal, doublecircuit, etc.). It will also be shown that this condition can be achieved with a singlereal transformation matrix T which is basically constant over the entire frequencyrange.2.2 Time Domain and Frequency Domain Solutions tothe MTL Equations2.2.1 FDLINE Solution to the MTL EquationsIn the FDLINE model, modal transformation matrices are used to decouple the N-coupled transmission lines to N-decoupled circuits. Each circuit is representedby a modal characteristic impedance Zcm and a modal propagation function γmcalculated the at each frequency fromZcm =√ZmYm, γm =√ZmYm (2.24)21At each decoupled circuit, the voltage and current waves are calculated inV (x) =Vf ke−γmx+Vbke+γmx (2.25)I(x) = I f ke−γmx+ Ibke+γmx (2.26)The voltage and current waves are related by the characteristic impedance inV (x) =Vf ke−γmx+Vbke+γmx = Zcm(I f ke−γmx+ Ibke+γmx)(2.27)Combining (2.25) (2.26) as in [69], we obtain the forward perturbation waveV (x)+ZcmI(x) = (Vk+ZcmIk)e−γmx (2.28)Note that in FDLINE [69], (2.28) can only be written in terms of modal de-coupled modes. Figure 2.3 illustrates FDLINE equivalent circuit at each mode. In(2.28), with x = 0 for the sending-end k of the line of Fig. 2.3, and x = ` forthe receiving end n, one obtains the voltage and current relationship for the lineequivalent circuit at node n,Vn+ZcmIn = (Vk+ZcmIk)e−γm` (2.29)The circuit of Fig. 2.3-right results from (2.29) after reversing the direction ofthe current. A similar analysis can be followed to obtain the equivalent circuit forthe line as seen from node k (Fig. 2.3-left).Ik (t) In (t)Ekh (t) Enh (t)k nZcm ZcmVk (t) Vn (t)+_ +_+_+_Figure 2.3: FDLINE modal equivalent circuit.The history sources in Fig. 2.3 are given byEnh = (Vk+ZcmIk)e−γm` , Ekh = (Vn+ZcmIn)e−γm` (2.30)22In FDLINE model, modal frequency dependent Zcm and e−γm` in the circuit ofFig. 2.3 are synthesized by means of BAF rational function approximation usingnegative poles and zeros to guarantee a passive realization.The BAF algorithm is based on an adaptation of the simple concept of asymp-totic fitting of the magnitude function, first introduced by Bode. The magnitude isplotted as a function of the logarithm of the frequency. The basic principle of theprocedure is to trace the approximated curve by straight-line segments. These seg-ments may be either horizontal or have a slope that is a multiple of 20 dB/dec. Thepoints of change of slope (corners or breaking points) define the poles and zeros ofthe rational approximating function [68].A major advantage of BAF over the other curve fitting algorithms which usecomplex poles (e.g., VF) is the process used to determine the time delay associatedwith the propagation modes [1] [69]. Because BAF uses only negative poles andzeroes, the resulting rational approximation is minimum phase-shift [92]. Thismeans that the real and imaginary parts are locked in a single function and the phaseis uniquely defined by the fit of the magnitude. For the propagation functions, therational fit has the same magnitude as the data function and the phase has a shiftof ωτ due to the time delay τ of the propagating mode. In FDLINE of [1] [69], thistime delay is determined by shifting the fitted function a for best correlation withthe original function. Since the original function is physical, it has to be minimum-phase-shift with a time delay and the BAF procedure assures the physical propertiesof the fit. This correlation procedure maintaining the physical constraints is notpossible when the fitting is not restricted to minimum-phase functions, like in thecase of VF [55].After approximating Zc by some poles and zeros, then it is simulated by a se-ries of RC parallel circuits. Using the trapezoidal integration rule, the RC networksare expressed in the form an equivalent resistance in series with a history voltagesource. Synthesis of the e−γ` results in two uncoupled parallel history voltagesource derived from the Inverse Laplace Transform of the partial fraction expres-sion of the approximated polynomials. The model can then be easily incorporatedinto a time domain solution simply in the form of a constant resistance in parallelwith a current source evaluated at each time step of the solution.Finally, the voltage and current of interest will be transferred from the modal23domain back to the phase domain using (2.21). As addressed in the literature,FDLINE approximates the transformation matrices in one frequency (e.g., 1 kHz)and takes only the real part of the matrices.The total operations count for FDLINE to simulate an N conductor system withPZc number of poles for fitting Zc and PAp number of poles for fitting e−γ` is givenbyN f dline =(N2+N∑1PZc)+(N2+N∑1PAp)(2.31)2.2.2 ULM Solution to the MTL EquationsThe ULM is a direct phase domain frequency dependent line model which is fun-damentally based on the idempotent decomposition (introduced as “column-wise”technique in [74]).Idempotent coefficients are matrices derived from the column-row multipli-cation of the modal transformation matrix and its inverse. For a function A, therelationship between the phase domain and modal domain are defined inAm = T−1AphT (2.32)Aph = TAmT−1 (2.33)Equation (2.33) can be rewritten in terms of the constituent columns for T andconstituent rows for T−1 for an N-phase multiconductor transmission lines inAph =[C1 C2 · · · CN]a11 0 · · · 00 a22 · · · 0....... . ....0 0 · · · aNNR1R2...RN=N∑i=1aii (CiRi) (2.34)The ULM uses VF algorithm to synthesize the modal line functions Yc and e−γ`and the idempotent coefficients matrices. VF method is a robust reformulationof the Sanathanan–Koerner iteration [103] using rational basis functions (partialfractions) instead of polynomials and pole relocation instead of weighting [57]. VF24is based on iteratively relocating an initial pole set to better locations [46]. Unlikethe BAF, VF requires both the magnitude and phase angle of the mode to be knownbefore the fitting can be done. Thus, the time constant for back-winding must bepre-calculated.In this model, VF was initially used to fit Yc column-wise using an identicalset of poles. Since Yc has no time delay, instead of fitting the sum of modes inphase domain using suitable sets of poles, fitting is done on the sum of the diag-onal elements of Yc. By taking the modal traveling time delays into account, thepropagation matrix is fitted in the phase domain. The resulting poles and time de-lays are then used for fitting e−γ` in the phase domain, under the assumption that allpoles contribute to all elements of e−γ`. The unknown residues are then calculatedby solving an over-determined system of linear equation as a least squares prob-lem. As all elements in e−γ` get identical poles, a column-wise realization [88]can be used which gives increased computational efficiency for the time domainsimulation.The total operations count for ULM to simulate an N conductor system with PYcnumber of common poles for fitting Yc and PAp number of poles for fitting e−γ` forN idempotents is given byNULM =(N2×PYc)+(N2×N∑1PAp)(2.35)2.2.3 NLT Solution to the MTL EquationsThe NLT implemented in this chapter follows the formulations in [75]. The NLTmethod solves the system directly in the frequency domain using nodal analysis.The current sources are converted from time to frequency, the admittance matrixis sampled in frequency, the system is solved, and finally the network voltages aretransferred from frequency to time. To relate the time and frequency domains, NLTrequires to define a time window of width Tc and a frequency window of width fcwith discrete samples. The time step size and the frequency step size are given by∆t =1fc, ∆ f =1Tc(2.36)The number of sample points (solution points) Ns is determined by25Ns =Tc∆t=fc∆ f(2.37)Two sampling scheme are introduced in [75] to discretize the parameters in thefrequency domain : regular sampling and odd sampling,∆ω = 2pik∆ f (regular) , ∆ω = pi (2k+1)∆ f (odd) (2.38)where, k = 0,1,2, . . . ,Ns−1.In addition to the discrete frequency, the Laplace variable s in the NLT includesa constant damping factor σ given in,s= σ + j∆ω (2.39)The factor σ is chosen to provide artificial damping to the time signal. Bydampening the time signal, the signal will fit into the time window Tc and the timealiasing problem due to the discretization of the admittance matrix Y(s) is reduced.There are three commonly used criteria for the calculation of the damping factor:Wilcox criteria σ1 [136], Wedepohl criteria σ2 [127], and error criteria σ3 [43]:σ1 = 2∆ω , σ2 =lnN2sTc, σ3 =− logεTc (2.40)where ε is the error tolerance.After identifying a time window and a frequency window, the input currentsources are transferred from time to frequency using the continuous Laplace trans-form.I(s) =∞∫0i(t)e−stdt∣∣s=σ+ j∆ω (2.41)Then, the system is solved using nodal analysis one frequency at a time inV(s) = Y−1(s)I(s) (2.42)In (2.42), V is the calculated vector of output node voltages, Y is the matrix ofadmittances evaluated at each frequency, and I is the injected current sources in theLaplace domain. To transfer the output network voltages from frequency to time,26NLT requires to use a boosting factor Cn and a weighting or windowing function(filter) W with the Inverse Fast Fourier Transform (IFFT) as shown inv(t) =ℜ(CnIFFT{V (s)W (s)}) (2.43)In (2.43), the real part of the signal is taken to eliminate the noise. The boostingfactor Cn is used to reverse the effect of the original artificial damping. Cn for theregular sampling and odd sampling is calculated inCn =1∆teσn∆t (regular) , Cn =2∆te(σ∆t+ jpiNs )n (odd) (2.44)where, n= 0,1,2, . . . ,Ns−1.In (2.43), windowing function W is usually applied to the NLT to attenuate theGibbs oscillations produced by the truncation of the frequency range. Lanczos andHanning are examples of such windowing functions [43].To incorporate multiconductor transmission lines in the NLT solution, a coupledexact-pi model are used as shown in Fig. 2.4.Y’ph2Z’phY’ph2Figure 2.4: Coupled exact-pi model. Frequency domain equivalent-circuit formulticonductor transmission line.To form this model, the series impedances Z′m and shunt admittances Y ′m ofthe decoupled exact-pi circuits are calculated in (2.45) by correcting the modal perunit length impedances and admittances (2.23) for the length ` and frequency. Thecorrection is performed by multiplying Zm and Ym by factors kz and ky in27Z′m = kzZm` , Y′m = kyYm` (2.45)kz =sinh(γm`)γm`, ky =tanh(γm`2)γm`2where γm is the modal propagation function calculated in (2.24).The coupled exact-pi model is obtained by transferring Z′m and Y ′m from modaldomain to phase domain inZ′ph = TVZ′mT−1I , Y′ph = TIY′mT−1V (2.46)It is important to note that the frequency variable to calculate the line parame-ters in the Laplace domain is s= σ + j∆ω .2.3 Case Studies and Simulation ConditionsIn this section, six line configurations are introduced to compare the solutions givenby FDLINE, ULM, and NLT:Test 1: Single-phase line.Test 2: Three-phase single-circuit horizontal line [66].Test 3: Three-phase single-circuit vertical line [34].Test 4: Three-phase single-circuit delta line [34].Test 5: Three-phase double-circuit one-tower delta line [66].Test 6: Three-phase double-circuit two-tower horizontal line [49].Single-circuit lines are considered untransposed to test FDLINE under asym-metrical line configurations and double circuit lines are chosen from the referencesthat previously criticized FDLINE for inaccuracy [66] [49]. In the test cases, onlyone conductor per bundle is considered in order to leave out the effect of the variousbundle modelling approaches used in different software packages.In Figs. 2.5 to 2.10 the line geometries and the terminal conditions of theconductors (open or shorted) are shown. The line length are selected based on thesource voltage level. As expected, the transients last longer in short lines due to thehigher perturbations caused by forward and backward waves travelling along the28short length.12 m1iscvoc1 Phase wires:d = 2.5 (cm),Rdc = 0.0576 (Ω/km),Ratio = 0.5,Source:Vs = 93.89cos(2π60t) (kV),Rs = 0.5 (Ω), Ls = 0.03 (H),Ground:ρ = 100 (Ω.m),µ = 1,Length: 150 (km),∆t = 5 µs.Figure 2.5: Single-phase transmission line (open and shorted).9.75 m17 m15.8 m1 2 3Phase wires: Bluebirdd = 4.475 (cm),Rdc = 0.0263 (Ω/km),Ratio = 0.364,Ground wires: 7No8 Awldd = 0.978 (cm),Rdc = 1.463 (Ω/km),Source: (@ each phase)Vs = 187.79cos(2π60t) (kV),Rs = 0.7 (Ω), Ls = 0.08 (H),Ground:ρ = 100 (Ω.m),µ = 1,Length: 240 (km),∆t = 40 µs.10 mi1v3Figure 2.6: Three-phase single-circuit horizontal transmission line.2918 m1 m7 m1237 m9 m Phase wires: Raild = 2.9591 (cm),Rdc = 0.0590 (Ω/km),Ratio = 0.375,Ground wires: 7No8 Awldd = 0.978 (cm),Rdc = 1.463(Ω/km),Source: (@ each phase)Vs = 281.69cos(2π60t) (kV),Rs = 1.2 (Ω), Ls = 0.13 (H),Ground:ρ = 100 (Ω.m),µ = 1,Length: 300 (km),∆t = 50 µs.i3v2Figure 2.7: Three-phase single-circuit vertical transmission line.23.3 m3.8 mPhase wires: Lapwingd = 3.82 (cm),Rdc = 0.03542 (Ω/km),Ratio = 0.375,Ground wires: 19No9 Awldd = 1.45 (cm),Rdc = 0.6821(Ω/km),Source: (@ each phase)Vs = 408.25cos(2π60t) (kV),Rs = 1.5 (Ω), Ls = 0.2 (H),Ground:ρ = 100 (Ω.m),µ = 1,Length: 480 (km),∆t = 80 µs.3.8 m5.7 m6 m125.2 m5.2 m35.5 mi3v1Figure 2.8: Three-phase single-circuit delta transmission line.3012.34 m7.01 mPhase wires: Grackled = 3.4 (cm),Rdc = 0.0472 (Ω/km),Ratio = 0.33,Ground wires: 7No5 Awldd = 1.39 (cm),Rdc = 0.7426 (Ω/km),Source: (@ each phase)Vs = 93.89cos(2π60t) (kV)Rs = 0.8 (Ω), Ls = 0.11 (H),Ground:ρ = 100 (Ω.m),µ = 1,Length: 100 (km),∆t = 16 µs.7.02 m8.17 m3.05 m8.23 m124514.33 m3 68.23 mi6v1Figure 2.9: Three-phase double-circuit one-tower delta transmission line.7.6 m15 m50 m9.0 m12 m13.6 m8.5 m2 34 5 6Phase wires:d = 3.2 (cm), Rdc = 0.05 (Ω/km),Ground wires:d = 1.0 (cm), Rdc = 0.5 (Ω/km),Ground: ρ = 100 (Ω.m), µ = 1,Length: 120 (km), ∆t = 20 µs.i1v2Vs = 93.89cos(2π60t) (kV),Rs = 0.5 (Ω), Ls = 0.03 (H),Vs = 187.79cos(2π60t) (kV),Rs = 0.7 (Ω), Ls = 0.08 (H),Figure 2.10: Three-phase double-circuit two-tower horizontal transmissionline.In the test cases, the transmission lines are connected to a balanced three-phasecosine source (Fig. 2.11) and the peak value of phase-a is applied at t = 0. Theequivalent source impedance corresponds to the impedance of the generator plusits step-up transformer. The conductors at the receiving-end of the line are eithershorted or open (connected to the ground with a resistance of 10−6Ω or 106Ω,31respectively). The simulations were run with the following conditions:456CoupledDouble-CircuitTransmission LineVsVst = 0123Rs LsRs LsTerminalsSingle-CircuitTransmission LineVst = 0 123Rs LsTerminals(a) (b)Figure 2.11: One-line diagram of the equivalent circuit for the tests. a)Single-circuit, b) Double-circuit.Simulation tools: Microtran v3.25 for FDLINE, PSCAD v4.5.2 for ULM, andMATLAB for NLT. Even though particular software packages were used to runthese tests, the cases are described with enough detail so that they can be run usingother EMTP software implementations.All methods: Shunt insulator losses Gins was taken as 3× 10−8 (S/km) [2].The simulations were run at the time steps ∆t’s indicated in the line configuration,Figs. 2.5 to 2.10. To guarantee the accuracy of the simulations, time steps are takento be 20 times smaller than the travelling time τ of the line with the length ` at thespeed of light υc = 3×108m/s(∆t = 120τ =`20υc).Line models: For the curve-fitting, maximum number of poles was set to 35over a frequency range from 10−2 to 107 Hz. The large frequency range is chosento increase the accuracy of fitting by allowing the asymptotic function to extendasymptotically beyond the range at which it will be used, which for switching tran-sients is about a few hundred kHz (for the calculation of line parameters, maximumfrequency is chosen at 1 MHz). For FDLINE, TV and TI were calculated at 1 kHz.NLT: Odd sampling, Error criterion (σ3) with ε taken as 10−3, and Hanningwindow. Frequency window fc is 1/∆t of the corresponding time steps. Numberof samples equals 220.322.4 Comparison of Simulation ResultsThe simulation results of the line energization tests, obtained with the NLT, FDLINE,and ULM solutions are shown in Figs. 2.12 to 2.17. Due to the closeness of the re-sults presented in each plot, the differences are best viewed using the glass magni-fier in a pdf viewer. The simulations are shown for a transient period t = 0 to t = 50ms. Steady-state errors are demonstrated in Table 2.1. Due to space limitations,only the phases with higher errors are shown.2.4.1 Single-Phase LineThe results for the single-phase line of Fig. 2.5 are shown in Fig. 2.12. Theseresults show that, in the absence of transformation matrices needed in the multi-phase cases, the BAF and VF curve-fitting algorithms used in FDLINE and ULM,respectively, can approximate the frequency dependence of the line functions quiteaccurately since the time domain solutions given by these line models perfectlymatch with the NLT solution which does not use these approximations.0 5 10 15 20 25 30 35 40 45 50-320-1600160320voc (kV)(a)NLT FDLINE ULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-1.4-0.700.71.4i sc (kA)(b)Figure 2.12: Simulation results for the single-phase line of Fig. 2.5.332.4.2 Three-Phase Single-Circuit LinesThe results for the single-circuit line tests are shown in Figs. 2.13 to 2.17. For thehorizontal line of Fig. 2.6, the open-circuit voltages at the receiving-end of the lineare shown in Fig. 2.13a. These results show that the solutions given by FDLINEand ULM are perfectly matched, and they follow the NLT solution very well (exceptin the jagged points of the curves in which the NLT deviates slightly). Also for theshort-circuit currents, Fig. 2.13b shows that the results of FDLINE and ULM for thehorizontal line are quite in agreement with the NLT solution.0 5 10 15 20 25 30 35 40 45 50-400-2000200400v3 (kV)(a)NLT FDLINE ULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-1.2-0.600.61.2i 1 (kA)(b)Figure 2.13: Simulation results for the horizontal line of Fig. 2.6.For the vertical line of Fig. 2.7, the open-circuit voltages at the receiving-endof the line are shown in Fig. 2.14a. These results show that FDLINE follows theNLT solution closer than ULM. The deviation of the open-circuit voltages at thepeak points and the slight phase shift drift of ULM can be more clearly observed.This phase shift is due to the lack of accuracy in ULM to determine the time delaywhich is caused by allowing non-minimum-phase shift fitting. For the short-circuit34currents, Fig. 2.14b shows that ULM is closer to the NLT solution than FDLINE. Thepeak value of short-circuit current of FDLINE is about 5% below the peak value ofthe NLT and ULM solutions. The reason for this difference might be due to takingthe transformation matrices at 1 kHz in FDLINE which for this simulation is not agood choice. This issue will be analyzed in Section 3.8.0 5 10 15 20 25 30 35 40 45 50-600-3000300600v2 (kV)(a)NLTFDLINEULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-2-101i 3 (kA)(b)Figure 2.14: Simulation results for the vertical line of Fig. 2.7.For the delta line of Fig. 2.8, the open-circuit voltages at the receiving-endof the line are shown in Fig. 2.15a. These results show that FDLINE and ULMfollow the NLT solution well. However, there is no perfect match between the threesolutions at the peak values of the few initial cycles. For the short-circuit currents,Fig. 2.15b shows that the short circuit currents of FDLINE and ULM match NLTsolution very well.2.4.3 Three-Phase Double-Circuit Line in the Same TowerThe FDLINE and ULM models performed very accurately for the single-phase lineand the three-phase single-circuit lines. As expected, however, larger differences350 5 10 15 20 25 30 35 40 45 50-1500-1000-50005001000v1 (kV)(a)NLT FDLINE ULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-2-1.5-1-0.500.51i 3 (kA)(b)Figure 2.15: Simulation results for the delta line of Fig. 2.8.between time-domain (FDLINE and ULM) and NLT simulations were observed forthe double circuit lines. The reason is due to more coupling between the conductorsin double circuit lines than the single circuit lines where the problem of frequencydependence of the transformation matrices for the line models and reaching theoptimal setting for the NLT solution become more challenging. The results of thesetests are shown in Figs. 2.16 to 2.17.For the double circuit delta lines of Fig. 2.9, the open-circuit voltages at thereceiving-end of the line are shown in Fig. 2.16a. As observed, the FDLINE andULM solutions are very close to each other and follow the NLT solution quite wellin the initial few cycles. However, the NLT solution deviates from the line modelsin the time span between 20 to 50 ms. For the short-circuit currents, Fig. 2.16bshows that both FDLINE and ULM solutions match the NLT solution perfectly.360 5 10 15 20 25 30 35 40 45 50-180-90090180270v1 (kV)(a)NLT FDLINE ULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-1.5-1-0.500.5i 6 (kA)(b)Figure 2.16: Simulation results for the double circuit delta line of Fig. 2.9.2.4.4 Three-Phase Double-Circuit Line in Separate TowersFigure 2.17 presents the results obtained for the two parallel horizontal lines ofFig. 2.10 mounted in two separate towers. For the open-circuit voltages, Fig. 2.17ashows that the solutions given by FDLINE and ULM are closer to the NLT solutionthan they were for the double circuit in the same tower case. This fact signifiesthat the choice of damping factor and window function for the NLT method weresuitable for this test. As observed, FDLINE is slightly closer to the NLT solutionthan ULM, particularly at the peak points. For the short-circuit currents, Fig. 2.17bshows that the results are very similar for both FDLINE and ULM. However, theNLT solution presents a slight vertical shift with respect to the results obtained withline models. As opposed to the open-circuit voltage, for the calculation of short-circuit current in two parallel horizontal lines, the NLT requires to use a differentdamping factor and/or windowing function so that it can match the result with theline models.370 5 10 15 20 25 30 35 40 45 50-320-1600160320v2 (kV)(a)NLT FDLINE ULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-2-1012i 1 (kA)(b)Figure 2.17: Simulation results for the two parallel horizontal line of Fig.2.10.2.4.5 Steady-state SolutionsTable 2.1 shows the steady state solutions obtained with NLT, FDLINE, and ULM inthe different test cases after the initial transient settles.In this table, the maximum error for the steady state voltages was 1.16% forFDLINE and 1.25% for ULM. The maximum error for the steady state currents was1.96% for FDLINE and 2.59% for ULM. Considering the complexity of frequencydependent line modelling, these are very good results for both FDLINE and ULM.2.5 General ObservationsOverall, simulation results showed that the FDLINE model gave similar results tothe ULM. These results indicated that, contrary to traditional belief, a constanttransformation matrix model like FDLINE is capable of representing multi-circuitasymmetrical line configurations.38Table 2.1: Comparison of steady state values of different solutions.Test kV/kA NLT FDLINE ULM1voc 167.81 167.9 0.05% 167.9 0.05%isc 1.155 1.159 0.36% 1.163 0.69%2v3 247.7 246.9 0.32% 248.0 0.12%i1 0.912 0.916 0.44% 0.919 0.77%3v2 397.5 396.9 0.15% 399.4 0.48%i3 0.970 0.951 1.95% 0.976 0.62%4v1 631.5 638.8 1.16% 633.6 0.33%i3 1.078 1.071 0.65% 1.072 0.56%5v1 100.6 101.1 0.50% 100.5 0.10%i6 0.809 0.803 0.74% 0.820 1.36%6v2 200.4 202.2 0.90% 197.9 1.25%i1 1.738 1.772 1.96% 1.783 2.59%The discrepancy between the time domain and frequency domain solutionsemerged in the simulation of double circuit lines where the solutions given byFDLINE and ULM were very close. The reason for the deviation of NLT from theline models might be due to inappropriate choice of error criterion for the dampingfactor and Hanning for the windowing function used for the NLT method. How-ever, this setting worked quite well for the single-phase line and three-phase single-circuit lines.Since the optimal set up for the NLT method is not easy to reach for non-expertusers, as a result, in the next chapter we will investigate a simpler methodology inthe frequency domain solutions which does not require the complicated set up asrequired for the NLT; while it provides very accurate results.39Chapter 3Discrete Time Fourier SeriesAlgorithmThis chapter derives a Discrete Time Fourier Series (DTFS) formulation for fre-quency domain solutions of fast electromagnetic transients in power systems. TheDTFS establishes a one-to-one mapping between Ns samples in a finite time win-dow and Ns samples in a finite frequency window. The system response can beevaluated using nodal analysis to obtain Ns/2 solutions in the frequency domainand this response is mapped back into Ns solution points the time domain. Theproposed methodology is based on the correct specification of the time windowand frequency window widths. The DTFS is simpler to implement than classicalNLT, which has been traditionally used for frequency domain solutions of powersystem transients. The DTFS can achieve similar levels of accuracy as the NLTwithout requiring user intervention to specify the value of the damping factor andfiltering windows required in the NLT.This chapter is organized as follows. First, the relationship between the dis-crete time domain and the discrete frequency domain is reviewed, then the stepsto correctly set up a DTFS circuit solution are presented. Next, guidelines are pro-vided to calculate a suitable time window and frequency window widths for theapplication of the DTFS. Simulation examples consisting of lumped R, L, C pa-rameters and frequency dependent transmission lines are introduced in preparationto compare the proposed DTFS with the conventional NLT and a reference EMTP40solution. Finally, the proposed DTFS is used to assess the accuracy of the FDLINEand ULM line models under asymmetrical line energizations and fault conditionsfor a double circuit vertical line.3.1 Relationship Between Discrete Time Domain andDiscrete Frequency DomainThe relationship between discrete time domain and discrete frequency domain in-troduced in this section follows the point of view of [67] where a mapping is estab-lished between time and frequency points. For a time window of Ns points there isa corresponding frequency window of exactly Ns points.Figure 3.1 illustrates a number of points in a time window of width Tc andits corresponding frequency window of width fc. As will be discussed in Section3.2.1, the width of the time window will be chosen based on the time constants ofthe circuit ([123],[140]). The width of the frequency window will be chosen basedon the maximum frequency that we want to simulate in the transient [28]. OnceTc and fc have been chosen, the number of sample points (solution points) Ns isdetermined byNs = Tc fc (3.1)The time step size and the frequency step size are given by∆t =1fc, ∆ f =1Tc(3.2)Equations (3.1) and (3.2) indicate that to simulate a larger time window for atime step ∆t we will need more solution points. This is exactly the same as in thetraditional EMTP criteria where one chooses the ∆t for a given maximum frequencyand then solves the system at Ns points (Ns = tmax/∆t) until tmax is reached.It should be noted that because the Fourier series to represent a real signal needsto be complex conjugates (Fig. 3.1), the maximum frequency that can be capturedin the solution is fNy (Nyquist frequency) which is one half of the sampling fre-quency,fNy =12∆t=fc2(3.3)41f∆xtfTcfcDTFSIDTFS-fNy= -5tfTcfc0+fNy= 5-1 1 2-2 3-3-4 4 0-1 1 2-2 3-3-4 4-fNy= -4.5 +fNy= 4.50 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8t∆Ns = 10 Ns = 9x(a) (b)|A| |A|0.5|A-5| 0.5|A+5|t∆f∆Figure 3.1: Correspondence of time domain and frequency domain for a dis-crete signal. a) Ns = even. b) Ns = odd.This means that there will be only Ns/2 distinct frequencies in the interval fromzero to+ fNy. The rest of the frequency components from− fNy to zero are complexconjugate of the corresponding positive interval.When simulating a physical system response, response samples are being pro-duced as the source applies input samples at t = n∆t, for n = 0,1,2, . . . ,Ns− 1.After the response to all these samples has been produced, the simulation timewindow will contain a snapshot that is mapped into a frequency window with fre-quencies in increments f = k∆ f , for k=−(Ns−1)/2, . . . ,−1,0,1, . . . ,(Ns−1)/2.Notice in Fig. 3.1a, for Ns=even that (Ns− 1)/2 is not an integer, then theclosest next integer is used. In Fig. 3.1b, for Ns=odd, the frequency component atfNy is split into two equal parts, one half on the negative side and the other half onthe positive side. Even though the frequency window of Fig. 3.1a shows 11 points,there are only 10 distinct frequencies since the − fNy and + fNy points are 1/2 ofthe same point and the coefficient at fNy is real (similarly, the coefficient at f = 0is also real). The frequency snapshot will provide frequency complex exponentialsthat can synthesize the values in the time snapshot. This synthesis is indicated in(3.4).The width of the time window Tc is determined by the desired simulation time42and how long the system response lasts. Similarly to the case of the frequency win-dow where only half of the window contains information, not all the time windowcontains information related to the solution of the transient. This is a crucial pointin the solution proposed in this research.The mapping between time and frequency in the DTFS is established as follows.A discrete-time signal x[n], having an arbitrary shape and a finite number of equallyspaced samples Ns, can be expressed as a sum of complex frequency componentssuch that when evaluated at t = n∆t we get the original x[n] values,x[n] =dNs−12 e∑k=b−Ns−12 cAkejk 2piNs n = IDTFS{Ak} (3.4)Equation (3.4) is the Inverse Discrete Time Fourier Transform (IDTFS) of thecoefficients of the Fourier series. These coefficients are the spectral coefficientsand are obtained fromAk =1NsNs−1∑n=0x[n]e− jk2piNsn = DTFS{x[n]} (3.5)Since the complex exponentials in the Fourier series are periodical functions,the synthesis of the time series is valid only for the original time span Tc.3.2 Power System Transient Solution Using the DTFSThe impulse response of the system h[n], which is the response of the system whena discrete-time impulse (a sample of value one) i[0] = 1 is applied, is illustrated inFig. 3.2a. If the sample of value one is applied at n = 1 instead of n = 0 (delayedimpulse), the shape of the response will be the same but shifted by one, i.e., h[n−1].For the next unit sample we will get h[n−2], etc. All these subsequent responseswill have the same shape as Fig. 3.2a but shifted to the right according to when thepulse is applied.If the magnitude of the input samples is i[n] instead of one, then the corre-sponding responses will be i[0]h[n], i[1]h[n−1], i[2]h[n−2], etc. The total outputby the time we reach point Ns will be the sum of all responses to previous inputsplus the response to the current input,43ht0 1 2 ... N1-1t1 2 3 ... N2-1zero completiontsimi30TcNs-1tset(a)(b)10.5...tsetFigure 3.2: Time window selection. a) System impulse response h[n], b) Cur-rent source i[n] (Dashed line: unit step with discontinuity, Solid line:unit step with the averaging of points at the discontinuity).v[n] = i[n]∗h[n] =Ns∑n=0i[n]h[Ns−n] (3.6)Equation (3.6) corresponds to the operation of convolution of the input functioni[n] and the system response h[n] evaluated at Ns. This procedure will allow us tofind the output v[n] knowing the input i[n] and the system response h[n].Suppose now that instead of obtaining the system response h[n] in the time do-main, our system is defined by a matrix of admittancesY [k] that has different valuesat different frequencies f = k∆ f . In this case, if we express the input source i[n] interms of its frequency components using (3.5) I[k] = DTFS{i[n]}, we can calculateV [k] by multiplying Y−1[k] by I[k] using the property that convolution in time (3.6)is equivalent to multiplication in frequency. OnceV [k] has been calculated for eachk, we can map the V [k] sequence to its corresponding v[n] sequence using (3.4).For the process described to work, we need a time window of adequate widthand a frequency window of adequate width so that the mapping between time do-main and frequency domain can be unique and can include all time dynamics andall frequency components in the simulated system.443.2.1 Defining the Time Window Width (Tc)Suppose we want to calculate the output v[n] to some input i[n] from n = 0 ton = N2− 1. Assume the system response h[n] (Fig. 3.2a) lasts for N1 samples.Then the time window Tc must be able to accommodate all the points of the source(N2) plus the length of the response h[n] (N1) to the last sample applied at N2−1.This gives us a length of N1+N2 which must fit in the time window Tc. Calling thelength of the source tsim and the length of the system impulse response tset , we willhave that Tc = tsim+ tset . This is the minimum width for the response modes of theimpulse response to die out within the time window. In Section 3.3.1, guidelinesare provided to calculate the correct tset .When transferring sequences between time window and frequency window, alltime sequences must have the same number of points in order to map the simulationNs-point time window into the simulation Ns-point frequency window. Since thesource sequence is only of length N2 we need to complete this sequence with zeros(“zero completion”) for the rest of the Tc window so that the system response cansettle without further source inputs. Since we stop the source at tsim, the simulationwill only be valid up to tsim. Tc can be made larger than the minimum required forthe purpose of increasing the frequency resolution (decreasing the spacing betweenfrequencies) without changing fc and at the expense of calculating more simulationpoints Ns. Adding extra zeros to increase the frequency resolution is called “zeropadding”. This is not to be confused with completing the source function withzeros to allow the system response to settle (“zero completion”).3.2.2 Defining the Frequency Window Width ( fc)For a given Tc we can increase the time resolution by decreasing the size of ∆t atthe expense of increasing the number of samples Ns = Tc/∆t. Decreasing ∆t willincrease the Nyquist frequency fNy and therefore higher frequency componentswill be simulated. As opposed to the criticality of the correct choice of the timewindow Tc, to incorporate all the transient and obtain correct results, the choice of∆t is determined by the maximum frequency we want to simulate in the transient.In discrete-time discrete-frequency analysis, if a transient contains frequencies upto 10 kHz and we choose a ∆t for fNy = 5 kHz, the system will be solved at discrete45frequencies up to 5 kHz. If the actual transient were to have higher frequencies,our solution will be only an approximation of the correct solution.In the case of using a continuous Laplace or Fourier transformation on a contin-uous source signal, higher source frequencies are captured. However, when solvingthe system for only a finite range of frequencies (of necessity because we can onlysolve for a finite number of points Ns), the time solution will present Gibbs oscil-lations due to the missing frequencies that we leave out of the solution.Guidelines for the calculation of fNy will be presented in Section 3.3.2.3.2.3 System Solution Using the DTFSAfter identifying a suitable time window where the transient solution will be lo-cated, and a suitable frequency window for the frequencies of interest, we cantransfer the source sequence from time to frequency using (3.5) and use nodal anal-ysis to solve the system one frequency at a time ( f = k∆ f ) for k= 0,1,2, . . . ,Ns/2,V[k] = Y−1[k]I[k] = H[k]I[k] (3.7)In (3.7), V is the calculated vector of output node voltages, Y is the matrix ofadmittances evaluated at each frequency (H is the inverse of Y), and I is componentk of the DTFS of the injected current sources in the window Tc. It is important tonote that the DTFS is obtained for each individual current source of the matrix i[n].To avoid division by zero when taking the inverse of Y, the DC frequency is takenas 10−4 Hz.After solving (3.7), the complementary part of the frequency spectrum of V [k](from − fNy to 0) is constructed by taking the complex conjugate of the positivehalf. For the NLT solution, this property only exists when regular sampling isconsidered (it does not apply to odd sampling). Finally, the time domain voltage atany node of interest is obtained by the synthesis of the time domain function fromits frequency components using the IDTFS of (3.4). After evaluating the IDTFS toobtain v[n] in the time window Tc, only the part of v[n] from zero to tsim representsthe output results. The rest of the time window corresponds to the de-energizationresponse and are “auxiliary points” that should be discarded.A source of errors when discretizing the input sources is related to disconti-46nuities in these sources. Consider, for example, a unit step function applied to asystem (Fig. 3.2b, dashed line). Because the coefficients in the spectral synthesisof the DTFS are continuous periodical functions, they cannot match discontinuitieswell. In addition, the DTFS formula uses a rectangular (forward Euler) approxima-tion of the “area” under the discrete time points, which results in a less accuratesynthesis than using, for example, a trapezoidal approximation to this area. Fig.3.16 in Section 3.5.2 illustrates the consequences of these problems. A consider-able improvement to this situation is to use the average of the values at the discon-tinuities, as indicated in Fig. 3.2b (solid line). This makes the calculated area equalto using the trapezoidal rule for the area under the discrete time points.3.3 Guidelines to Specify Tc and fcIn this section, guidelines are provided to calculate tset and fNy which are neededto calculate the width of the time and frequency windows, Tc and fc. The first stepis to extract the poles from the system’s transfer function. The n poles for a stablesystem of order n are:pi = αi± jβi , i= 1,2, . . . ,n (3.8)where, α and β correspond to the real part and the imaginary part of the systempole in this chapter.3.3.1 Calculation of the Time Window Width (Tc)As discussed earlier, the natural system responses must fit within the time window,and for this to happen Tc minimum should be tsim+tset . tsim is the simulation time ofinterest specified by user and tset is the time interval needed by the system responseto decay to zero. To ensure that the transient is completely damped within tset , werecommend to take this time as seven times the time constant of the slowest modeof the system τm (tset = 7τm). This choice will be more justified by performinga sensitivity analysis in Section 3.6. The negative of the real part of the poles isthe inverse of the time constants. The largest time constant of the system τm iscalculated asτm =1min(−αi) (3.9)47An approach to calculate τm for a system that includes frequency dependenttransmission lines will be presented in Section 3.4.1.3.3.2 Calculation of the Frequency Window Width ( fc)For an accurate representation of the system, the width of the frequency windowfc must be such that the system resonances fall within fNy = fc/2. As a guideline,fNy should be at least as large as the maximum of the following bandwidths:Br: bandwidth of the real part of the system poles,Bi: bandwidth of the imaginary part of the system poles,B`: bandwidth of the transmission line,Bs: bandwidth of the source,Bw: bandwidth of the switching for power electronic devices,Be: bandwidth of the estimation algorithm.The real part of the system poles corresponds to decaying exponentials in thetime domain. The decay of the exponentials was taken into account for the widthof the time window in (3.9). However, it also provides an index for the frequencybandwidth,Br = max(−αi)PPC (3.10)In (3.10), Br is augmented by a factor Points Per Cycle (PPC) to provide asuitable resolution. We recommend a PPC of 10.For some systems, the real part amplitude of some of the poles (α) might bemuch larger than that of other poles. Very large α means that the associated tran-sient decays very fast and for the time span of the simulation their contribution maybe negligible and ignoring them will prevent unnecessarily high values of fc. It isrecommended that values of α larger than 5×104 do not be taken into account tolimit the frequency window to 1 MHz and the corresponding ∆t to 1µs.The bandwidth corresponding to the imaginary part of the system poles aregiven byBi =max(|βi|)2piPPC (3.11)48For a system that includes transmission line, Br and Bi are calculated excludingthe line from the circuit. The bandwidth for a transmission line of length `, usingthe speed of light υc = 3×108m/s is calculated asB` =υc`PPC (3.12)Common sources in power systems are: sinusoidal, decaying exponential, andstep. For a sinusoidal source with frequency f0, Bs is calculated asBs = f0PPC (3.13)The bandwidth for the decaying exponential source is calculated using (3.10).The bandwidth for a step function of duration tsim that contains 99% of the energyis calculated from Parseval’s theorem [99],Bs =11tsim(3.14)For a system with power electronic devices such as inverters, Bw correspondsto the frequency of switching of power electronic switches.The estimation algorithms are usually used for power quality studies to extractthe amplitude and phase angle of different harmonics of a single [115] [117] [3].The bandwidth associated with these algorithms are taken as the highest order ofsignal’s harmonics.Since the DTFS provides a one-to-one mapping between the chosen time andfrequency windows, frequencies outside the specified frequency width do not exist.As a result, no frequency filtering windows are needed when using the DTFS, whilethey are normally needed when using the NLT to “chop” frequencies outside thiswindow. Also since the entire transient response of the system fits inside the chosentime window, the proposed DTFS method does not require, like the NLT, the use ofa damping factor σ in the system solution.3.4 Case StudiesThe test cases in this chapter include single-phase circuits and MTL-circuits.493.4.1 Single-phase CircuitsA 300 km single-phase transmission line is considered as a test case, as shown inFig. 3.3. A source is applied to the sending-end while the receiving-end is shortedthrough a resistance of 1 Ω. The conductor type is Rail and is placed 18 metersabove lossy ground with the resistivity of 100 Ω.m [116].18 mPhase wires: Raild = 2.9591 (cm)Rdc = 0.059 (Ω/km)Ratio = 0.375Ground: ρ = 100 (Ω.m), µ = 10=t-+Rs LsVsiscRsckm 300=lFigure 3.3: Example single-phase overhead transmission line.Fig. 3.4 shows three line models that can be used to represent the transmissionline of Fig. 3.3 in the frequency domain. The test circuits will consist of a voltagesource, a source impedance (Rs and Ls), one of the line models of Fig. 3.4, and aterminating short circuit resistance (Rsc). These circuits will be called: a) exact-picircuit, b) nominal-pi circuit, and c) RL circuit.(a) (b) (c)Rc Lc Rc LcC2C2Y’2Y’2Z’Figure 3.4: Single-phase line models. a) exact-pi , b) nominal-pi , c) RL.The purpose for introducing the circuits with lumped elements (Figs. 3.4b and3.4c) is for preliminary validation of the DTFS method with the reference analyticalformulas derived for these circuits.The line is energized at t = 0 via two types of voltage sources: a) a step voltagewith Vm = 1 kV for all three cases, and b) a cosine waveform with Vrms = 200 kVand 60 Hz for the exact-pi circuit. The equivalent source resistance and inductance50include the generator and step-up transformer and are taken as 1.2 Ω and 0.13 H,respectively.The exact-pi model of Fig. 3.4a is calculated at discrete frequencies (k∆ f fork = 0,1,2, . . . ,Ns/2) based on the procedure given in Section 2.2.3. The seriesimpedance and shunt admittances of the exact-pi circuit are obtained using the tra-ditional calculation for the per unit length line parameters presented in Section2.2.1. The shunt admittance Gins is considered as 2×10−9 S/km [39].Nominal-pi model (Fig. 3.4b) is a simplified version of the exact-pi model forshort lines and low frequencies. Under these conditions Gins is neglected, kz and kyare assumed to be one, and Z and Y are calculated at the power frequency. For theexample of Fig. 3.3, the lumped Rc, Lc, and C are calculated at 60 Hz as 35.4 Ω,0.68 H, 2.14 µF, respectively.The use of the nominal-pi model is introduced in the criteria to calculate theslowest time constant of the exact-pi model. The time constants of the exact-pimodel are mainly associated with the inductances and capacitances. Since induc-tances decrease with frequency and capacitances are assumed constant under TEMassumption, and the correction factors at very low frequencies are one, the slowesttime constant of the exact-pi model can be given by the nominal-pi model at verylow frequencies.The nominal-pi model can be further simplified into the RL model of Fig. 3.4cfor short lines by removing the shunt capacitances.The outputs are the short circuit current isc for all the three cases and the branchvoltage vL across the Lc in Fig. 3.4c. For comparisons, analytical formulas arederived for the output of the RL circuit and the nominal-pi circuit as shown inTable 3.1.3.4.2 MTL-CircuitsThe test case considered in this section is a double circuit vertical line with thegeometrical positions shown in Fig. 3.5 [34], which corresponds to the structurethat resulted in higher errors in [116].Modal lumped elements (Rc, Lc, C) of the corresponding nominal-pi model oftransmission line of Fig. 3.5 are given in Table 3.2. These parameters are calculated51Table 3.1: Analytical formulas for the RL and nominal-pi circuits.RL circuitisc(t) =VmReq(1− ept) , vL(t) = VmLcLeq eptReq = Rs+Rc+Rsc , Leq = Ls+Lc , p=−ReqLeq=−46.42Nominal-pi circuitisc(t) =−Vm4∑i=1kipi(1− epit)ki pi−1.2105×10−5 −9.3458×105+1.2345 −46.421−0.61723− j0.0061354 −8.1695+ j2926.3−0.61723+ j0.0061354 −8.1695− j2926.3at 0.01 Hz which are further used to calculate the slowest time constant of thecoupled exact-pi model.Table 3.2: Modal lumped elements of the nominal-pi circuit for the doublecircuit vertical line of Fig. 3.5 calculated at 10−4 Hz.Modes 1 2 3 4 5 6Rc (Ω) 17.7 17.7 17.7 17.7 17.7 17.7Lc (H) 4.3784 0.4925 0.4465 0.3631 0.3893 0.3882C (µF) 2.0747 4.5431 4.6674 5.7464 5.3475 5.5479The double circuit vertical line of Fig. 3.5 is energized at t = 0 for the followingtests: unbalanced faults, induced voltage, and step response. Because of spaceconsiderations, outputs are considered as the voltages and currents with highesterrors. The simulations for the short circuit currents gave much smaller errors than5218 m7 m7 m4.5 m123456Phase wires: Raild = 2.9591 (cm),Rdc = 0.059 (Ω/km),Ratio = 0.375,Gins = 2×10-9(S/km)Ground:ρ = 100 (Ω.m),µ = 1,Length: 300 km3P2Figure 3.5: Double circuit vertical transmission line.for the voltages.Unbalanced faults test:The equivalent circuit for this test is shown in Fig. 3.6. At the sending-end,each circuit is connected to a balanced three-phase cosine source and the peak valueof phase-a is applied at t = 0. The conductors 1, 3, 4, and 5 at the receiving-endof the line are open (connected to the ground with a resistance of 106 Ω), whereasconductors 2 and 6 are shorted with a resistance of 1Ω. The outputs are the voltageat open conductor 3 and current at the shorted conductor 6.Induced voltage test:Figure 3.7 illustrates the system diagram for the induced voltage test. Thistest can be performed for the maintenance of individual three-phase circuits. Inthis example, the bottom circuit is grounded at both sending-end and receiving-endwith 10 Ω resistances. The top circuit is energized with a three-phase balancedcosine voltage at sending-end with open conductors 1 and 3 and shorted conductor2 with a 1 Ω resistance at the receiving-end. The outputs are the current at shortedconductor 2 and the voltage at the terminal of conductor 6.530.13 H456CoupledDouble-CircuitTransmission LineVsVst = 0123+v3_0.13 H1.2 Ω1.2 ΩVs :Vs1 =Vs4 = 281.7cos(120πt)Vs2 =Vs5 = 281.7cos(120πt-2π/3)Vs3 =Vs6 = 281.7cos(120πt+2π/3)km 300=l1 Ω1 Ωi6Figure 3.6: System diagram for the unbalanced faults test.Vs :Vs1 = 281.7cos(120πt)Vs2 = 281.7cos(120πt-2π/3)Vs3 = 281.7cos(120πt+2π/3)123456CoupledDouble-CircuitTransmission Linekm 300=lt = 0i2+v6_0.13 H1.2 Ω1Ω10 Ω10 ΩVsFigure 3.7: System diagram for the induced voltage test.Step response test:The equivalent circuit for this experiment is shown in Fig. 3.8. In this test,conductor 1 is energized with a step voltage with the amplitude of 1 kV at sending-end, and is kept open at the receiving-end. The rest of the conductors are shortedwith a 10 Ω resistance at both sending-end and receiving-end. The voltage at openconductor 1 and the voltage at the shorted conductor 4 are considered as outputs inthis test.For simplicity, the test systems of Figs. 3.6 to 3.8 are called “MTL-circuits”.3.5 Simulation ResultsIn this section, the simulation results of the proposed DTFS method are comparedwith the NLT method, a reference EMTP solution run at a very small ∆t = 1 µs,54123456CoupledDouble-CircuitTransmission Linekm 300=lt = 0+v1_123+v4_10 Ω10 Ω10 Ω10 Ω0.13 H1.2 Ω1 kVFigure 3.8: System diagram for the step response test.and the analytical formulas for the nominal-pi and the RL circuits. The analyticalformula is used to assess the accuracy of the DTFS and NLT in a form of an absoluteerror. For the exact-pi circuit and MTL-circuits, the solutions given by the DTFSand NLT are benchmarked with two EMTP line models: the FDLINE model usingMicrotran v3.25 and the ULM using PSCAD v4.5.2. In the curve-fitting process ofFDLINE and ULM, the maximum number of poles is set to 35, the frequency rangeconsidered is from 10−2 to 107 Hz (as discussed in Section 2.3, the maximum rangefor frequency is chosen at 100 MHz to increase the accuracy of the curve-fitting).For FDLINE, the transformation matrices are taken at 1 kHz (for MTL-circuits). TheNLT method used for these comparisons follows the regular sampling formulationgiven in [75]. For benchmarking purposes, the same time and frequency windowsproposed in Section 3.3 for the DTFS are also used for the NLT. The dampingfactors used for the NLT are: Wilcox criteria σ1, Wedepohl criteria σ2, and errorcriteria σ3 with ε taken as 10−6. The windowing functions for the NLT solution areconsidered as Lanczos, Hanning, and Tukey with 2Q/N = 0.3 [11].The poles of the nominal-pi circuit at very low frequencies (10−4 Hz) are usedto calculate the largest time constant of the exact-pi circuit. These poles are:−4.6729× 105,−16.446,−5.0496± j2006.8. To calculate the largest time con-stant of the MTL-circuits, modal lumped elements indicated in Table 3.2 are usedto form six decoupled circuits for each line energization tests. Using the eigen anal-ysis, the smallest amplitude real part pole for MTL-circuits is calculated as 2.58202for mode 1-circuit for all of the three test systems.55Table 3.3: Time-frequency parameters for the DTFS.Parameters RL circuit Nominal-pi circuit Exact-pi circuit MTL-circuitτm (s) 0.021543 0.122406 0.198035 0.387294Tc (s) 0.200801 0.906842 1.436245 2.761058∆ f (Hz) 4.980055 1.102728 0.696260 0.362180Br (kHz) 0.464198 0.464210 0.169231 0.169231Bi (kHz) - 4.657351 - -B` (kHz) - - 10 10Bs (kHz) 0.22 0.220.22 (step) 0.22 (step)0.6 (cosine) 0.6 (cosine)fc (kHz) 0.928396 9.314702 20 20Ns 187 8447 28725 55222∆t (µs) 1073.80 107.36 50 50Table 3.3 shows the step by step implementation of the guidelines of Section3.3 to calculate the Tc and fc windows for the DTFS method.For all methods, the simulation time of interest is considered as 50 ms (un-less otherwise stated). τm is calculated from the knowledge of the system poles(poles of the RL and nominal-pi circuits are shown in Table 3.1). The time windowdetermined for the DTFS is tsim+ 7τm. The corresponding frequency resolution is∆ f = 1/Tc. For the frequency window, Br and Bi are calculated from the knowl-edge of the system poles, B` is determined by the length of the transmission line,and Bs depends on the shape of the waveform. For the nominal-pi circuit, we calcu-late Br by taking the second largest amplitude of the real part of the pole. (We skipthe largest real part amplitude pole based on the threshold introduced in Section3.3 to limit the smallest time step to 1µs.) The Nyquist frequency fNy is chosenas the largest of the above bandwidths. This gives fc = 2 fNy. The total number ofsamples in the time domain and frequency domain is Ns = Tc fc. The correspondingtime step size is ∆t = 1/ fc.For DTFS and NLT with regular sampling the system is solved from 0 to fNy(k= 0 to Ns/2). The negative frequencies half of the spectrum is built by taking the56complex conjugate of the positive frequencies half. As discussed, the simulationresults are valid from t = 0 to t = tsim for the DTFS and, for comparisons purposes,only this interval will also be considered for the NLT. The results for a number ofcase studies are presented next.3.5.1 RL Circuit SimulationFigure 3.9 shows the short circuit current isc calculated with all methods for thestep voltage source for a time window Tc = 50 ms.0 5 10 15 20 25 30 35 40 45 50Time (ms)0510152025i sc (A)DTFSNLT <1NLT <2NLT <3EMTPAnalytical formulaFigure 3.9: isc for the RL circuit energized with a step voltage.As observed in Fig. 3.9, the DTFS solution matches very well the NLT, EMTP,and analytical formula. For the simulation time of interest tsim= 50 ms, the solutiongiven by the NLT did not present any numerical oscillations at the tail of the signal,as a result, the NLT does not require to use windowing function for this simulation.Maximum absolute error for the DTFS is 0.26 A and for the NLT (for the threechoices of damping factor presented here) is 0.13 A. This error occurs in the firstsample for both solutions which is due to the use of forward Euler integrationmethod in the conventional IDTFS. Ignoring the first few samples of the simulation,the average absolute error for the DTFS is 0.006 A and for the NLT is 0.001 A. Thisfact signifies that the NLT solution can be very accurate when the optimal choice ofdamping factor is considered, but only as long as the correct time window width is57also taken into consideration.This simulation was run for tsim = 0.5 s. Figure 3.10 shows the close up ofthe comparison of DTFS with NLT and EMTP at the tail of the signal. It can beobserved that the NLT solution contains oscillations for tsim = 0.5 s, while for thesame simulation time the solution given by the DTFS follows the analytical formulaand the EMTP solution very well. Therefore, NLT (unlike the DTFS) requires to usethe windowing function to reduce the errors.0.4 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5Time (s)2626.226.426.626.827i sc (A)DTFSNLT <1NLT <2NLT <3EMTPAnalytical formulaFigure 3.10: Close up of isc for the RL circuit energized with the step voltagefor tsim = 0.5 s.Figure 3.11 shows the effect of different windowing functions to reduce theerrors of NLT solution in Fig. 3.10. The plot shows the result for NLT with σ3(similar results were observed when σ1 and σ2 were selected).It can be observed that Lanczos’ window caused minor attenuation to the os-cillations since they still show up after the filtering. The filtering effect of Hanningwindow is stronger than the Lanczos’ window, however, it leaves a DC offset afterthis process which impairs the accuracy of the NLT solution. Tukey window has abetter filtering impact with minor superimposed oscillations around the referenceanalytical formula. However, the DTFS provides a more accurate solution than theNLT without the need to use any windowing functions for this simulation.Figure 3.12, shows the results for the branch voltage vL in the RL circuit for astep voltage source for t from zero to tsim = 50 ms with the same set up considered580.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5Time (s)26.5826.5926.626.6126.62i sc (A)DTFSNLT <3 + LanczosNLT <3 + HanningNLT <3 + TukeyEMTPAnalytical formulaFigure 3.11: Effect of different windowing functions to reduce the oscilla-tions of the NLT solution in Fig. 3.10.in Fig. 3.11. vL(kV)Figure 3.12: vL for the RL circuit energized with the step voltage.It can be observed that when no filter is applied to the NLT, the solution presentssome superimposed Gibbs oscillations, more noticeably, in the beginning and tailportions of the simulation (magnified), while the results for the DTFS, EMTP, andanalytical formula are perfectly matched with each other and do not present theseoscillations. The cause of error in the NLT is the chopping of frequencies that59results in Gibbs phenomena. These results illustrate the need for a windowingfunction in the NLT, while, the proposed DTFS formulation does not require thisfiltering. The maximum absolute errors for the NLT solution with the applicationof different windowing functions are: 38 V for Tukey, 42 V for Lanczos, and 83V for Hanning, while the maximum absolute error for the DTFS is 0.7 V whichis significantly smaller than those for the NLT. Applying Hanning’s window tothe NLT in this case resulted in an absolute error larger than when no windowingfunction was used (81 V). This shows the sensitivity of NLT to the windowingfunction used.Ignoring the first few samples, the average absolute error for the NLT is 0.3 Vfor Tukey, 0.66 V for Lanczos, and 0.98 V for Hanning, while this value for DTFSis 0.2 V. The error analysis performed for the simulation shows that for the NLT togain higher accuracy, we need to search for other choice of damping factors andwindowing functions than the conventional ones used in this study.The scenario of Fig. 3.12 was also used to illustrate the effect of odd samplingin the transient response of NLT. Figure 3.13 compares NLT with odd samplingwith DTFS, EMTP, and analytical formula for tsim = 0.5 s. vL(kV)Figure 3.13: Test of NLT with odd sampling in the simulation of vL for theRL circuitAs opposed to DTFS and NLT with regular sampling which solved the systemup to Ns/2, using the odd sampling for the NLT requires to solve the system for60the total number of samples Ns. This happens, as discussed in Section 3.1, due tolosing the complex conjugate property in the frequency domain when odd samplingis considered. Therefore, the number of operations is doubled when NLT uses oddsampling.As observed in Fig. 3.13, the NLT solutions deviate from the EMTP and ana-lytical formula at some point of the simulation while the DTFS perfectly coincideswith the reference solutions until tsim is reached. Similar to the Gibbs oscillations,these deviations for the NLT depend on the value of the damping factor. In thissimulation, σ1 by causing less error offers a better setting for the NLT. Figure 3.13also shows that NLT underestimates the transient peak by 42 V for all the dampingfactors.Figure 3.14 illustrates the impact of different windowing functions to eliminatethe deviation of NLT method in Fig. 3.13. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5Time (s)00.20.40.60.81vL(kV)DTFSNLT 1 + LanczosNLT 1 + HanningNLT 1 + TukeyEMTPAnalytical formulaFigure 3.14: The effect of filtering on the accuracy of NLT with odd sampling.Simulation of vL for the RL circuit.This effect shows up as an overestimation, unlike the results of Fig. 3.13, ofthe transient peak by 34 V for Tukey, 5.6 V for Hanning, and 4.2 V for Lanczos,while the maximum error of DTFS is 0.7 V.613.5.2 Nominal-pi Circuit SimulationFigure 3.15 shows the simulation results for isc for the nominal-pi circuit energiza-tion with a step voltage source and tsim= 50 ms. The results show that the solutionsgiven by all methods are basically identical to each other. The maximum absoluteerrors are: 4 mA for DTFS, 1 mA for NLT with σ1, and 0.1 mA for NLT with σ2 orσ3, which are all negligible. The NLT did not require to use windowing functionfor this simulation.0 5 10 15 20 25 30 35 40 45 50Time (ms)0510152025i sc (A)DTFSNLT <1NLT <2NLT <3EMTPAnalytical formulaFigure 3.15: isc for the nominal-pi circuit energized with the step voltage.The scenario of Fig. 3.15 was also used to test the DTFS algorithm for the stepfunction input with discontinuities (dashed line in Fig. 3.2b). The solution givenby the DTFS without the proposed averaging method is denoted as “DTFS0”. Fig.3.16 is the magnification of Fig. 3.15 for the simulation span between 40 to 45 msin which DTFS0 is compared to the rest of solutions. As observed in Fig. 3.16,DTFS0 drifts slightly from the rest of the solutions.3.5.3 Exact-pi Circuit SimulationFigures 3.17 and 3.18 show isc for the exact-pi circuit calculated with the DTFS,NLT, FDLINE, and ULM for the simulation time of 50 ms. Figure 3.17 shows iscfor a step voltage source, whereas Fig. 3.18 shows the result for a cosine sourceenergization. The results of Figs. 3.17 and 3.18 confirm that the DTFS and NLT6240 40.5 41 41.5 42 42.5 43 43.5 44 44.5 45Time (ms)22.222.422.622.82323.223.423.6i sc (A)DTFSDTFS0NLT <2EMTPAnalytical formulaFigure 3.16: Details of Fig. 3.15 illustrating the effect of the first-last pointcorrection.coincide well with the FDLINE and ULM time domain models in the EMTP.0 5 10 15 20 25 30 35 40 45 50Time (ms)05101520253035i sc (A)DTFSNLTFDLINEULMFigure 3.17: isc for the exact-pi circuit energized with the step voltage.An incorrect selection of the time window can lead to aliasing problems in theNLT and to accuracy problems for the proposed DTFS method. Even though theissue of aliasing does not apply to the DTFS because it is a one-to-one mappingbetween frequency and time, if the system response does not have time to dieout within the time window, the DTFS will not include the response to the last630 5 10 15 20 25 30 35 40 45 50Time (ms)-1.5-1-0.500.511.5i sc (kA)DTFS NLT FDLINE ULMFigure 3.18: isc for the exact-pi circuit energized with the cosine voltage.samples of the source. The example of the exact-pi circuit energization by a stepvoltage source with an incorrectly chosen time window can be used to illustrate thisproblem. Suppose the RL circuit is used to calculate the time constant for the exact-pi circuit at very low frequencies instead of the proposed nominal-pi circuit. The RLcircuit gives a time constant τRL = 0.06 s, which is smaller than the time constantobtained with the nominal-pi circuit (τNpi = 0.198 s). Based on the proposed criteriafor the calculation of Tc in Section 3.3.1, choosing smaller time constant results inreducing the width of the time window.Figure 3.19 shows a magnification of Fig. 3.17 in which the solution given bythe DTFS-τRL (DTFS with incorrect Tc) is compared to the DTFS-τNpi (DTFS withcorrect Tc), NLT, FDLINE and ULM. Choosing an incorrect time window results inan inaccurate results for the DTFS-τRL, while the results for DTFS-τNpi is perfectlymatched with the rest of solutions.3.5.4 MTL-Circuits SimulationFor the simulation of MTL-circuits, difficulties were experienced to set the timewindow Tc, damping factor, and sampling scheme for the NLT method.As mentioned earlier, for benchmarking purposes, the same Tc for the DTFS wasconsidered for the NLT solution. The Tc consists of tsim= 0.05 s and tset = 7τ = 2.71640 0.5 1 1.5 2 2.5 3Time (ms)01234i sc (A)DTFS-=N:DTFS-=RLNLTFDLINEULMFigure 3.19: Error in the DTFS due to incorrect selection of the time window.s which results in 2.76 s for MTL-circuits. For the same size of Tc, the NLT solutioncould match with the DTFS, FDLINE, and ULM solutions only when σ1 for bothregular/odd sampling or σ3 for odd sampling were selected. Since the operationof regular sampling with the use of complex conjugate property is faster than theodd sampling (89.2 s compared to 171.3 s for the simulation of Fig. 3.6), σ1 withregular sampling was considered as the set up for NLT. For other settings to work,the NLT method required to consider a larger time window, as for example: 3.2 sfor σ2 and 2.79 s for σ3 for regular sampling and 3.12 s for σ2 for odd sampling.Windowing functions had minor impact to enhance the result of the NLT for itsoptimum setting with σ1 or to reduce the size of the typical time windows for othersettings.These observations indicate that obtaining accurate results with the NLT re-quires experience to choose a correct damping factor. While DTFS is much simplerto implement and does not require this setting.Figure 3.20 shows the voltage at open conductor 3 of the double circuit verticalline in the unbalanced faults test of Fig. 3.6. In this simulation, the results obtainedwith the DTFS is compared with NLT, FDLINE, and ULM for tsim= 50 ms. As can beobserved in Fig. 3.20, the solutions given by DTFS and NLT with σ1 are perfectlysuperimposed, and these solutions are closely followed by the FDLINE and ULM650 5 10 15 20 25 30 35 40 45 50Time (ms)-600-3000300600v3 (kV)DTFS NLT <1 NLT <2 FDLINE ULMFigure 3.20: v3 in the unbalanced faults test of Fig. 3.6.line models. However, the solution given by NLT with σ2 does not agree with therest of the matched results. This conflict shows the sensitivity of the NLT to thesize of the time window and the choice of the damping factor. For the rest of thesimulations, σ1 and regular sampling will be considered as the set up for the NLTmethod.In Fig. 3.21, short circuit current i6 for the double circuit vertical line in theunbalanced faults test of Fig. 3.6 is shown. These results indicate that the solutionsgiven by FDLINE and ULM line models are suitably in agreement with the DTFSand NLT frequency domain solutions.Figure 3.22 shows the simulations for v6 and i2 in the induced voltage test ofFig. 3.7 for tsim = 50 ms. As can be observed, for both v6 and i2 the DTFS isperfectly matched with the NLT, and both solutions are closely followed by theFDLINE and ULM line models.660 5 10 15 20 25 30 35 40 45 50Time (ms)-2.5-2-1.5-1-0.500.511.5i 6 (kA)DTFSNLTFDLINEULMFigure 3.21: i6 in the unbalanced faults test of Fig. 3.6.0 5 10 15 20 25 30 35 40 45 50-4.5-3-1.501.53v6 (kV)(a)DTFS NLT FDLINE ULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-1.201.22.4i 2 (kA)(b)Figure 3.22: v6 and i2 in the induced voltage test of Fig. 3.7.67Figure 3.23 shows the solutions obtained with DTFS, NLT, FDLINE, and ULMfor the open conductor 1 for the double circuit vertical line in the step response testof Fig. 3.8. In this figure, the results of the DTFS and NLT are suitably superim-posed. However, time domain models have minor differences with respect to thefrequency domain solutions at the peaks. The maximum error for FDLINE is 17 Vand for ULM is 47 V.0 5 10 15 20 25 30 35 40 45 50Time (ms)-0.500.511.522.5v1 (kV)DTFS NLT FDLINE ULMFigure 3.23: v1 in the step response test of Fig. 3.8.Simulation of the step response test of Fig. 3.8 is also used to illustrate theconsequence of not correcting the transformation matrices for the modal switching.Figure 3.24a shows the real part of TI(ω) for the double circuit vertical line ofFig. 3.5 which contains the modal switching for some elements at 400 Hz and 2.5kHz. In Fig. 3.24b, TI(ω) was corrected for the modal switchings with the use ofNewton-Raphson algorithm [135].The DTFS was used as an example in this section to address the modal switch-ing problem. However, this is not a typical problem of the DTFS and can also occurfor other time and frequency domain solutions. The DTFS solution which is notcorrected for the modal switching is denoted as “DTFS-MS”.Figure 3.25 shows the simulations for the voltage at shorted conductor 4 ofdouble circuit vertical line in the step response test of Fig. 3.8. Results are obtainedwith DTFS-MS (contains modal switching), DTFS (corrected for modal switching),NLT, FDLINE, and ULM.6810-2 10-1 100 101 102 103 104 105 106-0.8-0.400.40.8Real of TI(!)(a)10-2 10-1 100 101 102 103 104 105 106Frequency (Hz)-0.8-0.400.40.8Real of TI(!)(b)Figure 3.24: The real part of TI(ω) for the double circuit vertical line of Fig.3.5. a) with modal switching. b) corrected for modal switching.0 5 10 15 20 25 30 35 40 45 50Time (ms)-1-0.500.51v4 (V)DTFS DTFS-MS NLT FDLINE ULMFigure 3.25: v4 in the step response test of Fig. 3.8 illustrating the effect ofmodal switching problem.69As observed in Fig. 3.25, modal switching appears as high frequency oscil-lations in the time domain result, more noticeably at the beginning of the signal.Similar to the previous results, the DTFS and NLT solutions are best coincided.FDLINE follows the frequency domain solutions closer than ULM in this simula-tion. The maximum absolute error for FDLINE is 1.0154 kV and for ULM is 1.0841kV compared to the DTFS as a reference.Overall, by agreeing with the DTFS and NLT solutions, these tests also showthat the frequency dependent line models FDLINE and ULM make a very good job infitting the frequency dependent wave functions in developing their correspondingtime domain models.3.6 Sensitivity of the DTFS to the Width of the tsetIn this section, a sensitivity analysis is performed to show the accuracy and com-putational cost of the DTFS method based on the guidelines presented in Section3.3.1. Figure 3.26a illustrates a relative error of the DTFS method with respect tothe reference solution for the width of tset ranging between 0.1 to 10 s. In Fig.3.26b, the computational efficiency of different simulations are illustrated. Thetest cases were selected from the examples in Section 3.4 which resulted in highererrors:– isc in the RL circuit,– isc in the nominal-pi circuit,– isc in the exact-pi circuit energized with the step voltage,– v3 in the double circuit vertical line in the unbalanced faults test.The reference solutions for the error analysis are the analytical formulas for theRL and nominal-pi circuits and the solution given by the EMTP line models withinthe first travelling time (which is zero as shown in Fig. 3.19). The errors are thennormalized by the steady state values of the analytical formulas for circuits withlumped elements and EMTP line models for the frequency dependent lines. In Fig.3.26a, arrows indicate the width tset = 7τm as proposed in Section 3.3.1.As can be observed in Fig. 3.26a, DTFS has large error when tset is close to0.1 s. With the proposed set up for tset , the errors are as follows: 0.7% for the RLcircuit, 0.003% for the nominal-pi circuit, 0.006% for the the exact-pi circuit, and703.5× 10−6% for the double circuit line. Taking tset = 2τm for the double circuitline test results in 1.7×10−5% error.By comparing Figs. 3.26a to 3.26b we can conclude that increasing the width oftset from the proposed criterion will increase the computational cost at the expenseof no significant gain in the accuracy. 10-110010110-610-410-2100102Error (%)RL CircuitNominal- CircuitExact- CircuitDouble Circuit Line10-1100101tset(s)10-2100102104Computational Time (s)Figure 3.26: Sensitivity of the DTFS method to the width of the tset for differ-ent tests. a) Relative error, b) Computational time.As discussed in Section 1.1, frequency domain solutions are not as efficientas time domain solutions in terms of the computational cost. Their application islimited to verify the accuracy of the corresponding time domain models.713.7 General ObservationsIn this chapter, we proposed a simple methodology to implement a frequency do-main solution based on the DTFS for power system transient analysis. The proposedalgorithm specifies the time window width based on the slowest time constant ofthe system and the frequency window width based on the maximum frequency ofinterest in the transient.The accuracy of the proposed DTFS was verified by comparing the simulationresults with the analytical formulas obtained for circuits with lumped elements.The DTFS solution coincided well with the FDLINE and ULM in the simulation ofa double circuit vertical line under asymmetrical line energization and terminalconditions.Simulation results showed the sensitivity of the NLT solution to the selectionof damping factor, windowing function, and time window width. For some sim-ulations in this chapter, similar to the simulation results in Chapter 2, we couldnot obtain good results with the NLT that compare well with the DTFS and EMTPsolutions due to the difficulty to adjust the NLT parameters. While, the DTFS usesa simpler implementation than the NLT which facilitates its usage for non-expertusers. Sensitivity analysis performed in Section 3.6 verified the accuracy and com-putational efficiency of the proposed guidelines for the application of DTFS.Simulation results for the double circuit vertical line in Section 3.5.4 repeatedthe conclusion in Chapter 2 that FDLINE, by approximating the transformation ma-trices at one frequency, is very accurate for the simulation of asymmetrical lineconfigurations. In the next section, we further investigate the effect of calculatingthe transformation matrices at different frequencies on the accuracy of the FDLINE.3.8 Sensitivity of FDLINE to the Frequency ofApproximation for the Transformation MatricesThe suggested frequency for the calculation of the transformation matrices forFDLINE is 1 kHz [2]. In this section, we perform a sensitivity analysis to findthe effect of this frequency on the accuracy of FDLINE.The double circuit vertical line of Fig. 3.5 in the unbalanced faults test of Fig.3.6 is chosen as an example for this test. FDLINE is simulated with transformation72matrices calculated for a frequency range from 10−2 Hz to 106 Hz with 10 pointsper decade.The accuracy of the FDLINE is evaluated by calculating a relative error withrespect to the DTFS solution for the transient period and steady state condition.The transient error is calculated as the average of absolute error of FDLINE andDTFS over a period of 0.1 s where most transient happens. The transient error isnormalized by the steady state value of DTFS solution, and it is then scaled by afactor 100. The steady state error is calculated as a relative error of FDLINE andDTFS at tsim = 0.6 s where the transient is almost settled. The transient and steadystate errors of ULM are also calculated for comparison purposes. Since ULM doesnot approximate the transformation matrices, these errors are constant over thefrequency range of study. Figure 3.27 shows the transient and steady state errorsfor FDLINE and ULM with respect to the DTFS.As observed in Fig. 3.27, FDLINE has the highest errors for frequencies smallerthan 1 Hz. For the frequency range of 10 Hz to 10 kHz (including 1 kHz as thetraditional setting) FDLINE can be more accurate than ULM for all the simulations,except for i6 that ULM is slightly better (0.88% error for ULM versus 1.22% errorfor FDLINE).From the results of Fig. 3.27, one can conclude that choosing the frequency atabout 60 Hz can result in minimum transient and steady state errors for FDLINE.For simplicity, we refer to the solution obtained with the traditional setting at 1kHz as “FDLINE” and for the new setting at 60 Hz as “FDLINE*”.In Fig. 3.27, the maximum transient error for the open voltages are 3.72% forULM, 3.00% for FDLINE, 2.69% for FDLINE*, and for short circuits are 2.13% forFDLINE, 0.94% for ULM, 1.22% for FDLINE*. The maximum steady state errorsfor both open voltages and short circuit currents are 1.31% for ULM, 0.81% forFDLINE, and 0.71% for FDLINE*. These results indicate that FDLINE can be moreaccurate when the frequency of approximation for the transformation matrices arecalculated at 60 Hz than the traditional setting at 1 kHz. These results were alsoobserved for a variety of other line configurations.7310-2 100 102 104 106 10801020Error of V1 (%) Transient ErrorFDLINEULM10-2 100 102 104 106012Steady-state Error10-2 100 102 104 106 108024Error of i 2 (%)10-2 100 102 104 106024610-2 100 102 104 106 10801020Error of V3 (%)10-2 100 102 104 10602410-2 100 102 104 106 10801020Error of V4 (%)10-2 100 102 104 106012310-2 100 102 104 106 108051015Error of V5 (%)10-2 100 102 104 106024610-2 100 102 104 106 108Frequency (Hz)024Error of i 6 (%)10-2 100 102 104 106Frequency (Hz)0246Figure 3.27: Sensitivity of FDLINE to the frequency of approximationSimulation results of Figs. 3.20–3.23 are magnified in Figs. 3.28–3.30 to illus-trate the improvement of FDLINE with the new setting.74Figure 3.28 shows the close up for v3 and i6 for the unbalanced faults test ofFig. 3.6. Results clearly show that FDLINE* follows the DTFS solution closer thanFDLINE and ULM for both open voltage v3 and short circuit i6.8 9 10 11 12 13 14 15 16Time (ms)100180260340420v3 (kV)16 17 18 19 20 21 22 23 24 25Time (ms)-400-300-200-1000100v3 (kV)5.5 6 6.5 7 7.5 8 8.5Time (ms)-2.3-2.2-2.1-2-1.9i 6 (kA)DTFSFDLINE: 1 kHzFDLINE: 60 HzULM14 14.5 15 15.5 16 16.5 17 17.5Time (ms)0.50.540.580.62i 6 (kA)Figure 3.28: Close up of v3 and i6 in the unbalanced faults test of Fig. 3.6.75Details of i2 and v6 for the induced voltage test of Fig. 3.7 are illustrated inFig. 3.29. As observed in this figure, FDLINE* is slightly closer than FDLINE, andboth models are closer than ULM when compared to the DTFS solution.7.5 8 8.5 9 9.5 10 10.5 11Time (ms)1.861.891.921.95i 2 (kA)17 17.5 18 18.5 19Time (ms)-0.72-0.7-0.68-0.66i 2 (kA)5 7 9 11 13Time (ms)-4.3-3.8-3.3-2.8-2.3v6 (kV)15 16 17 18 19 20Time (ms)0.51.11.72.3v6 (kV)DTFSFDLINE: 1 kHzFDLINE: 60 HzULMFigure 3.29: Close up of i2 and v6 in the induced voltage test of Fig. 3.7.76Figure 3.30 shows the transient solutions for v1 in the step response test of Fig.3.8 for two periods: a) 0< t < 20 ms and b) 0.1< t < 0.15 s. As can be observed,for the simulation time of interest up to 20 ms, the solutions obtained with FDLINEand ULM are basically matched with the DTFS solution. In this period, FDLINE isnot sensitive to the frequency of approximation. However, as simulation advances,FDLINE* continues to follow the DTFS with minimum errors in the peaks, whileULM and FDLINE solutions include a noticeable phase shift drift with larger errorsat the peaks. The deviation of the ULM is larger than FDLINE which is mainly dueto sacrificing the minimum-phase shift condition in the fitting functions.0 4 8 12 16 20Time (ms)-0.250.250.751.251.752.25v1 (kV)0.1 0.11 0.12 0.13 0.14 0.15Time (s)0.850.951.051.15v1 (kV)DTFS FDLINE: 1 kHz FDLINE: 60 Hz ULMFigure 3.30: v1 in the step response test of Fig. 3.8.77Chapter 4Revised MulticonductorTransmission Line Equations,Proposed Frequency-DependentLine ModelThis chapter proposes additional physical constraints on the classical MTL equa-tions. The main conceptual premise is that the voltage wave and the current waveshould be collocated and travel with the same propagation function along the line.This leads to the analysis of the way the resistances of the ground return path areincluded in the series impedances. A set of revised equations is proposed RMTLthat has the same propagation function for currents and voltages. As opposed tothe MTL equations that require complex frequency dependent transformation ma-trices for their diagonalization, the RMTL equations can be diagonalized using asingle real constant transformation matrix. Based on the RMTL equations, a new“frequency dependence simulation” line model “FDLM” is proposed that is muchsimpler than ULM, while still being fully frequency dependent. The validity ofFDLM is verified with comparisons with FDLINE, ULM, and a recently developedDTFS frequency domain solution.784.1 Inconsistencies in the Classical MTL EquationsAs discussed in Chapter 2, the per unit length series impedance and shunt admit-tance matrices are traditionally calculated inZ = (Rc+Re)+ jωLloop (4.1)Y = Gins+ jωCextThe wave propagation constants γV and γI are calculated from the series impedanceand shunt admittance of (4.1),γV =√ZY =√(Rc+Re+ jωLloop)(Gins+ jωCext) (4.2)γI =√YZ =√(Gins+ jωCext)(Rc+Re+ jωLloop)(4.3)In general, the product of two matrices is not commutative, and it is assumedin the classical MTL solution (e.g., [126][93]) thatZY 6= YZ (4.4)If (4.4) were to be true, thenγV 6= γI (4.5)In a conceptual experiment, let us assume that we apply an excitation to amulticonductor uniform semi-infinite transmission line. For this line, there willbe no return waves and the voltage and current propagation equations (2.15) and(2.16) will only have forward components,V f (x) = V f ke−γVx , I f (x) = I f ke−γIx (4.6)where V f (x) and I f (x) are related as in (4.7).V f (x) = ZcI f (x) (4.7)A fundamental assumption of TEM propagation, on which the MTL equationsare based, is that in a long uniform transmission line, there are no waves propagat-79ing perpendicular to the line and the longitudinal propagation equations (4.8) aresufficient to describe the line.∂ 2V∂x2= (ZY)V ,∂ 2I∂x2= (YZ)I (4.8)This assumption (e.g., [93]), requires that skin effect in the capacitances beignored and that C remains constant at its electrostatic value Cext . This is theassumption made in the existing EMTP line models.The ratio, in a matrix sense, between the voltage wave and the current wave isthe characteristic impedance matrix Zc of (4.9).Zc = Y−1√YZ (4.9)Since the line is assumed to be uniform and infinite, the Zc that relates thevoltage and current waves at any point x of the line should be the same. However,this is not enforced in the classical MTL equations. Relating in (4.6), V f (x) andI f (x), with the assumption that γV 6= γI , it follows thatV f (x) =[Zce−(γV−γI)x]I f (x) (4.10)In this case then the relationship between the forward voltage wave and theforward current wave becomes a function of the position x, which contradicts thephysical assumption of a semi-infinite uniform line, where we should be able toapply a given voltage at any point x and get the same current. Propagation accord-ing to (4.10) makes the voltage and current waves separate from each other as thewave advances, that is, to lose “collocation”. The problem of enforcing, at leastthe property of symmetry in the product of the characteristic admittance and thepropagation function is studied in [108]. In this study, however, we go beyondsymmetry, and require that the characteristic impedance does not change anywherein the semi-infinite line and, therefore, that γV − γI in (4.10) be zero.4.2 Collocation of Current and Voltage WavesTo resolve the physical inconsistency of (4.10), we postulate that for a uniformsemi-infinite transmission line in an isotropic medium, (4.2) and (4.3) must be80equal to each other, regardless of the mathematical argument that the product oftwo matrices is in general (but not necessarily) non-commutative. If we impose thecondition γV = γI , we obtain, for the forward wave in (4.10), and with a symmetri-cal argument for the corresponding backward wave,V f (x) = ZcI f (x) , Vb(x) =−ZcIb(x) (4.11)In (4.11) current and voltage are always related by the same Zc, independentlyof where the waves are in the line. We refer to (4.11) as the wave collocationcondition. For this condition to be true, it is required that the matrix product of(4.4) be commutative, that is,ZY = YZ (4.12)Next, we will look into factors that can be used to enforce the condition (4.12)in the MTL equations.4.3 Effect of the Earth Return in the Matrix ofImpedancesAs discussed in Chapter 2 and with reference to Fig. 2.2, self and mutual elementsof the complex loop inductance L′loop are calculated inL′loop11 =µ2piln2(h1+ p)r1eq, L′loop12 =µ2pilnD12′d12(4.13)The imaginary part of the complex loop inductances of (4.13) are used to cal-culate the ground return resistances,Re11 = jω(jℑ(L′loop11))=−ωℑ(L′loop11) (4.14)Re12 = jω(jℑ(L′loop12))=−ωℑ(L′loop12)If the height of conductor 2 in Fig. 2.2 is different from the height of conductor1, the values for self and mutual loop inductances would be different in (4.13). Asa resultRe11 6= Re12 (4.15)81Other methods (e.g., [13], [118]) can be used to calculate the values of the Re’s,however, the values obtained are also different for the selves and the mutuals.Notice that Re11 and Re12 in (4.14) have to be physically the same, even thoughthe formulas will give slightly different results.To illustrate this conflict, we consider the test circuit of Fig. 4.1 for the line ofFig. 2.1 which shows two conductors 1 and 2 and ground return (it is assumed thatthe other conductors are open). Assume ∆x = 1 m to determine elements Z11 andZ12 in Ω/m for the series impedance matrix.I1x∆x yEV1xgx gyReRc Lloop11I1Rc Lloop12012V2xV1yV2y+_Figure 4.1: Test circuit to determine the self and mutual impedances in theseries impedance matrix Z.Assume that at position y both conductors are shorted to ground and that atposition x, a source E is connected to conductor 1, while conductor 2 is open, Inthis test circuit, current I1 flows in conductor 1 and returns through ground. Withrespect to the voltage drop equations (2.1), we will haveV1x−V1y = E = Z11I1 , V2x−V2y =V2x = Z12I1 (4.16)and from the circuit solution,Z11 = (Rc+Re)+ jωLloop11 (4.17)Z12 = Re+ jωLloop12 (4.18)Since in Fig. 4.1, conductor 2 is not carrying any current, the only resistancecontributing to Z12 is the resistance of the earth circuit Re encountered by I1. No-82tice also that in this experiment, Re in (4.18) is the same Re in (4.17) because thecorresponding voltage drop is due to the same current I1 flowing in the same earthpath.We can repeat the tests in Fig. 4.1 for the other elements of the impedancematrix. Assume, for example, that we swap the terminal conditions in conduc-tors 1 and 2, and we inject a current into conductor 2 and measure the voltage inconductor 1. In this experiment, we will get Z22 and Z21.If we now require, as in the first experiment, that the Re for the mutual Z12 bethe same as for Z22, this will make Z12 different from the value Z21 obtained inthe first experiment. This will contradict the condition of symmetry in an isotropicmedium, that is, Z12 = Z21.Here, to be consistent with the description on N conductors with a commonground return in the MTL equations, we will make all values of Re equal to theaverage of all contributions to Z at each frequency Re and still use the slightlynon-physical approximations of (4.13) for the calculations.Figure 4.2 shows elements of the earth resistance matrix Re for the doublecircuit vertical line of Fig. 3.5 for a frequency range 10−2 Hz to 106 Hz. In Fig.4.2, Re given by (4.14) for the MTL equations are illustrated with green curves,while the black curve is the average of the green curves to be used by the RMTLequations. The maximum difference between the elements of Re in MTL model areused to calculate an error which is shown with the red curve in Fig. 4.2. This erroris normalized by the average value of all contribution of Re and further expressedin percentage. As can observed in Fig. 4.2, for the normal range of transients inpower systems 1 to 10 kHz, these differences are within 10-25%.By factoring out the identical frequency dependent Re from the matrix of groundresistance, the remaining matrix can be represented as “all-ones” matrix for an N-conductor system [1]N×N . With this description, matrix of impedance Z in theRMTL model can be expressed asZ = Rc+Re[1]+ jωLloop (4.19)where Rc is a diagonal matrix, [1] is an N×N matrix of ones, and Lloop is afull matrix.8310-2 10-1 100 101 102 103 104 105 106Frequency (Hz)10-610-410-2100102104Re (+ /km) in MTLRe (+ /km) in RMTLError (%)Figure 4.2: Comparison of Re in MTL and RMTL models for the double cir-cuit vertical line of Fig. 3.5.4.4 Validity of the Collocation Condition ZY = YZThe plots in Fig. 4.3 show the validity of the wave collocation condition for thedouble circuit line of Fig. 3.5 thatZY(YZ)−1 = I (4.20)where, I refers to the identity matrix in this chapter. Ideally, the diagonalsshould be one and the off-diagonals should be zero.The curves in Fig. 4.3 are obtained using the conventional Z and Y parametersfor the MTL equations and also with the ground resistances Re taken as equal forall entries in the Z matrix for the RMTL equations.Both sets of curves are very close, except that enforcing the collocation con-dition improves (reduces) the imaginary part of all elements of the matrix (Fig.4.3-bottom). As it can be seen in the plots, the maximum errors in the diagonalsare about 2% and in the off-diagonals about 5%.These results suggest that it is the nature of the system to satisfy the collocationproperty, and that there might be some physical inconsistencies in the equations orin the determination of the parameters of these equations that prevent the collo-cation condition to be exact. The new RMTL equations proposed in this chapter8410-2 10-1 100 101 102 103 104 105 1060.980.9911.011.02Real of ZY(YZ)-1 Diagonal MTLRMTL10-2 10-1 100 101 102 103 104 105 106-0.06-0.0300.030.06Real of ZY(YZ)-1 Off-Diagonal10-2 10-1 100 101 102 103 104 105 106Frequency (Hz)-0.03-0.02-0.0100.010.020.03Imaginary of ZY(YZ)-1All ElementsFigure 4.3: Validity of the collocation condition. If ZY = YZ, the prod-uct ZY(YZ)−1 should be one for the diagonals and zero for the off-diagonals.enforce the collocation condition as exact and incorporate some corrections to themodelling.4.5 Diagonalization of the RMTL EquationsAs indicated in (2.22), to diagonalize the MTL equations, it is necessary to findthe matrix of eigenvectors TV that diagonalizes the product ZY, and the matrix ofeigenvectors TI that diagonalize the product YZ. In general, since it is assumedthat the products ZY and YZ and are different and complex, TV and TI will alsobe complex and different from each other. Also, due to the frequency dependenceof the elements of Z, TV and TI will be frequency dependent. From the theory in85[56] [128], the TV and TI matrices that diagonalize ZY and YZ will also diago-nalize the individual Z and Y matrices, Zm = T−1V ZTI and Ym = T−1I YTV , wheresubscript m denotes the modal diagonalized quantities. The proposed RMTL equa-tions are based on satisfying the collocation condition (4.12) ZY=YZ. To achievethis, while minimizing the effect of the inconsistencies of the MTL description, wepropose that instead of first diagonalizing the matrix products ZY and YZ, as nor-mally done, we first diagonalize the individual Z and Y matrices. We will start byfinding a transformation matrix T that diagonalizes Z, and then we will use T todiagonalize Y with a very good approximation. The accuracy of this approxima-tion will be assessed in Section 4.6 and Section 4.7. This way the products ZY andYZ will be diagonalized by the same matrix, that isZphYph = TZmT−1TYmT−1 (4.21)= TZmYmT−1 = TYmZmT−1 = YphZphBecause T will diagonalize the matrix of impedances, we will additionally re-quire that T preserves the nature of the real and imaginary parts of Z, that is, thelosses and time constants should be the same in the original phase system and in themodal system. The same applies to T diagonalizing Y. This leads to the argumentthat there must be a solution for T that is real. This argument is consistent withthe property in matrix algebra that the eigenvectors that diagonalize a symmetri-cal positive definite real matrix have a real solution (in addition to other possiblecomplex solutions).4.6 Transformation Matrix T to Diagonalize Z and YSince Z is subjected to skin effect, while Y is formulated based on the constantGins and Cext , it is a more restricted problem to diagonalize Z than to diagonalizeY. We will then first try to find a transformation matrix that diagonalizes Z, with agood degree of accuracy and with certain desired properties, while keeping in mindthat this matrix must also diagonalize Y.Now, suppose in (4.19) that instead of diagonalizing the full Z directly, we firstdiagonalize its dominant component Lloop. Since Lloop is a symmetrical positivedefinite real matrix, there is a real solution for its diagonalizing matrix T. Thissolution is shown in Fig. 4.4 for the line of Fig. 3.5.8610-2 10-1 100 101 102 103 104 105 106-0.8-0.400.40.8Real of T10-2 10-1 100 101 102 103 104 105 106Frequency (Hz)-0.4-0.200.20.4Imaginary of TFigure 4.4: RMTL formulation. T matrix to diagonalize Lloop for the line ofFig. 3.5. The same matrix T diagonalizes ZY and YZ.In Fig. 4.4, we can see that T is basically constant up to about 1 kHz, withsome variations after this frequency. The maximum variance is about 2× 10−5for the frequency range between 10−2 Hz and 104 Hz, and slightly larger at 10−4for the range from 104 to 106 Hz. This is in contrast with the strong frequencydependence of the complex TV and TI matrices of the current MTL formulation.For reference, these matrices are shown in Fig. 4.5 calculated using MATLAB andusing the method of [135] to avoid modes switching.To proceed with the diagonalization of (4.19), we will take a single value of Treal. Here we use the value at 100 Hz which is asymptotic in the low frequenciesregion. Any other frequency in this region gives basically the same results.We now apply T−1(.)T to each term on the right hand side of (4.19). Since Rcis already diagonal, we will get the same diagonal Rc. Next, we apply T−1(.)T tothe matrix of all ones [1]. Notice that a matrix of all ones is a special case of a8710-2 10-1 100 101 102 103 104 105 106-0.8-0.400.40.8Real of T(!)10-2 10-1 100 101 102 103 104 105 106Frequency (Hz)-0.4-0.200.20.4Imaginary of T(!)TV(! ) TI(! )Figure 4.5: MTL formulation. Real and imaginary part of complex TV and TImatrices to diagonalize ZY and YZ for the line of Fig. 3.5.symmetrical matrix where all the self and mutual elements are one. Using sym-metrical components or other akin transformation gives, for exact diagonalization,Zsel f +(N− 1)Zmutual = N for the zero sequence, mode 1, and Zsel f −Zmutual = 0for the line modes, modes 2 to N. Next, the diagonalizing matrix T is applied tothe Y matrix in (4.1). Since Gins is already diagonal, the concern is to diagonalizethe matrix of capacitances Cext . Table 4.1 shows the errors in the diagonalizationof the matrix of ones [1] and Cext , for the double-circuit vertical line of Fig. 3.5.Notice that since T is constant, and the matrix of ones and Cext are also con-stant, the errors in Table 4.1 are the same for all frequencies. As seen in this table,the diagonalization is practically exact for modes 3, 4, and 6. For the matrix ofones, the errors are 0.04% for mode 1, 0.0011% for mode 2, and 0.23% for mode5. For Cext , the errors are 0.9% for mode 1, 0.3% for mode 2, and 0.09% for mode5. This accuracy has been verified for other line configurations with similar results.88Table 4.1: Errors in the diagonalization of the matrix of ones [1] and Cext forT calculated at 100 Hz.modes 1 2 3 4 5 6[1]Exact 6 0 0 0 0 0Approx. 5.9978 1.1×10−5 0 0 0.0023 0Cext (pF/m)Exact 3.4574 7.5722 7.7790 8.9125 9.2466 9.5774Approx. 3.4890 7.5491 7.7792 8.9123 9.2380 9.5773The results shown in Table 4.1 justify the validity of diagonalizing Z and Y sepa-rately with the same matrix T. Next we look at the capability of T to diagonalizethe products ZY and YZ.4.7 Transformation Matrix T to Diagonalize ZY and YZThe premise in obtaining a single real constant T was that this matrix would beable to diagonalize both ZY for the propagation of the voltages and YZ for thepropagation of the currents, and thus satisfy the collocation condition (4.12).In the plots of Fig. 4.6 we show the magnitude of the off-diagonal elementsafter diagonalization. Ideally, these elements should be zero. To show the relativeerrors, the values were divided by the geometrical average of the correspondingdiagonal elements at each frequency. The largest errors occur between 10 and 100Hz for the off-diagonal elements (5,6) and (2,6). However, since in the line modelthe off-diagonals are made exactly zero, these errors are not significant as long asthe values of the diagonals are accurate.In 4.2 and 4.3, we assess the accuracy of the diagonal elements resulting fromthe diagonalization of ZY using the proposed T matrix versus the traditional TVmatrix. For each frequency, the first row corresponds to T and the second rowto TV . The accuracy is similar for the diagonalization of the propagation of thecurrents (not shown). The numbers in these tables have been divided by the speedof light. Even though all numbers are negative, they are shown as positive in thetables to save space.Notice that the real parts (4.2) are one to two orders of magnitudes larger thanthe imaginary parts (4.3). The errors in the real parts are less than 0.1% for all89 (6,2) (5,6) Rest of the off-diagonal elements (6,5) (2,6) Figure 4.6: Diagonalization of T−1ZYT and T−1YZT. The plots show theabsolute value of the off-diagonal elements.modes, except at very low frequencies where the magnitudes are very small. Theerrors for the imaginary parts remain less than 1.5% for all modes for frequen-cies below 100 kHz (which is the normal range for switching transients in powersystems), while they grow up to 4% for the ground mode, and higher for some ofthe line modes, for higher frequencies. For frequencies above 1 MHz, the basicassumptions of TEM propagation become questionable [124].90Table 4.2: Real part of the eigenvalues of ZY calculated using the proposedsingle T versus exact diagonalization (all values are negative).f Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 61 Hz3.418e−7 1.026e−7 1.017e−7 1.016e−7 1.022e−7 1.019e−73.370e−7 1.067e−7 1.017e−7 1.019e−7 1.016e−7 1.038e−7100 Hz2.685e−3 1.378e−3 1.369e−3 1.367e−3 1.372e−3 1.370e−32.699e−3 1.367e−3 1.369e−3 1.370e−3 1.367e−3 1.368e−31 kHz2.155e−1 1.354e−1 1.345e−1 1.339e−1 1.342e−1 1.340e−12.168e−1 1.343e−1 1.346e−1 1.340e−1 1.339e−1 1.339e−110 Hz1.725e+1 1.338e+1 1.331e+1 1.324e+1 1.325e+1 1.324e+11.736e+1 1.328e+1 1.332e+1 1.324e+1 1.323e+1 1.323e+1100 Hz1.471e+3 1.327e+3 1.323e+3 1.319e+3 1.319e+3 1.319e+31.477e+3 1.321e+3 1.324e+3 1.318e+3 1.318e+3 1.318e+31 MHz1.367e+5 1.320e+5 1.319e+5 1.317e+5 1.317e+5 1.317e+51.369e+5 1.328e+5 1.319e+5 1.317e+5 1.317e+5 1.317e+5Table 4.3: Imaginary part of the eigenvalues of ZY calculated using the pro-posed single T versus exact diagonalization (all values are negative).f Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 61 Hz4.615e−7 8.454e−7 8.708e−7 9.961e−7 1.032e−6 1.070e−64.607e−7 8.458e−7 8.707e−7 9.961e−7 1.070e−6 1.033e−6100 Hz4.110e−4 9.195e−5 9.468e−5 1.083e−4 1.130e−4 1.163e−44.131e−4 9.169e−5 9.528e−5 1.079e−4 1.162e−4 1.117e−41 kHz3.331e−2 2.174e−3 2.235e−3 2.560e−3 2.720e−3 2.750e−33.348e−2 2.260e−3 2.325e−3 2.522e−3 2.730e−3 2.632e−310 Hz2.341e0 6.377e−2 6.540e−2 7.491e−2 8.245e−2 8.050e−22.377e0 7.760e−2 8.113e−2 7.443e−2 8.007e−2 7.713e−2100 Hz1.191e+2 1.974e0 2.019e0 2.313e0 2.643e0 2.485e01.231e+2 3.626e0 4.586e0 2.339e0 2.472e0 2.386e01 MHz4.590e+3 0.621e+2 0.634e+2 7.258e+1 8.471e+1 7.800e+14.806e+3 1.571e+2 2.304e+2 7.592e+1 7.772e+1 7.513e+14.8 Proposed Frequency-Dependent Line ModelThe wave collocation conditions (4.11) (4.12) result in simpler RMTL wave propa-gation equations. For convenience, we repeat (4.12) here,91ZY = YZ (4.22)With this condition, (2.15) (2.16) becomeV(x) = V f ke−γx+Vbke+γx (4.23)I(x) = I f ke−γx+ Ibke+γx (4.24)where the same propagation constant applies to both, voltages and currents,γ =√ZY =√YZ (4.25)The voltage and current waves are related byV(x) = V f ke−γx+Vbke+γx = Zc(I f ke−γx+ Ibke+γx)(4.26)where the characteristic impedance is given by (4.9) and repeated here,Zc = Y−1√YZ (4.27)Equations (4.22) to (4.27) constitute the RMTL and will be used to formulatethe proposed FDLM model.Combining (4.23) (4.24) as in [69], we obtain the forward perturbation waveV(x)+ZcI(x) = (Vk+ZcIk)e−γx (4.28)Note that in the existing FDLINE model of [69], (4.28) can only be written interms of individual decoupled modes (scalar quantities), while in FDLM, becauseof having a single γ , (4.28) applies also to the full phase-coordinates matrices.Equation (4.28) leads to the equivalent circuit of Fig. 4.7. With x = 0 for thesending end k of the line, and x= ` for the receiving end n, one can write,Vn+ZcIn = (Vk+ZcIk)e−γ` (4.29)The circuit of Fig. 4.7-right results from (4.29) reversing the direction of thecurrent. A similar analysis can be followed to obtain the circuit of Fig. 4.7-left for92Ik (t) In (t)Ekh (t) Enh (t)k nZc ZcVk (t) Vn (t)+_ +_+_+_Figure 4.7: FDLM equivalent circuit in the frequency domain. Quantities canbe matrices in the phase domain or scalar quantities in the modal do-main.the line as seen from node k. The history sources in Fig. 4.7 are given byEnh = (Vk+ZcIk)e−γ` , Ekh = (Vn+ZcIn)e−γ` (4.30)For multiconductor lines in phase coordinates, the quantities in Fig. 4.7 arematrices. For decoupled modal coordinates, the quantities are scalars for eachpropagation mode.The modal (diagonal) series impedance and shunt admittance matrices aregiven byZm = T−1ZphT , Ym = T−1YphT (4.31)where now in FDLM we have a single real frequency-independent transforma-tion matrix T in (4.31) instead of two complex frequency dependent transformationmatrices TV and TI in (2.23).The modal characteristic impedance Zcm and the modal propagation functionγm are calculated at each frequency fromZcm =√ZmYm, γm =√ZmYm (4.32)In FDLM, like in FDLINE, the modal frequency dependent characteristic impedanceand propagation functions are synthesized using negative poles and zeroes for apassive realization. The improved BAF of [1] is used for the synthesis. For the lineof Fig. 3.5, the maximum errors with BAF were less than 0.1% for the propagationfunctions and less than 0.5% for the characteristic impedance functions over theentire frequency range.934.9 Numerical Results: Transient SimulationsThis section compares transient simulation results obtained with the proposed FDLMline model and the traditional EMTP FDLINE and ULM formulations. In previouswork in [116], we showed that FDLINE and ULM gave very similar accuracy forasymmetrical overhead line configurations.The analysis in the present section further makes the case for the validity ofusing a single transformation matrix. Here, the proposed FDLM model is comparedto FDLINE and ULM for the double circuit vertical line of Fig. 3.5 [34], whichcorresponds to the structure that resulted in higher errors in [116].To validate the representation of the transformation matrices, a very accurateDTFS solution was implemented as proposed in Chapter 3.Simulation tools: PSCAD v4.5.2 for ULM, Microtran v3.25 for FDLINE, andMATLAB for FDLM and DTFS.All methods: simulation time of interest = 50 ms, time step ∆t = 50µs.DTFS: sampling frequency fc = 20 kHz, time window width Tc = 3 s.Line models: The number of poles for the fitting of the line functions was setto a maximum of 35 for FDLINE, ULM, and FDLM. An average of about 10 polesper function was used by the models. The frequency range was from 10−2 to 107Hz. For FDLINE, TV and TI were calculated at 1 kHz. For FDLM, T was calculatedat 100 Hz.Figures 4.8 and 4.9 show the results for the open circuit voltage v3 and the shortcircuit current i6 for the test of Fig. 3.6 using FDLM, FDLINE, and ULM for the timedomain solutions and DTFS for the frequency domain solutions. FDLINE and ULMuse the MTL equations while FDLM uses the RMTL equations. Two solutions arecalculated with DTFS: a) For the MTL equations, calculating TV and TI at eachfrequency, and b) For the RMTL equations, calculating T at each frequency.Figure 4.10 shows the instantaneous errors over the time line of the simulationuntil steady state is reached. In general the errors for the currents were smaller thanthe errors for the voltages.In Fig. 4.9 and Fig. 4.10 for v3, we can observe that all curves are basicallysuperimposed around the first (negative) peak. On the second (positive) peak, ULM940 5 10 15 20 25 30 35 40 45 50-600-3000300600v3 (kV)(a)0 5 10 15 20 25 30 35 40 45 50Time (ms)-2.4-1.201.2i 6 (kA)(b)DTFS: RMTLFDLMDTFS: MTLFDLINEULMFigure 4.8: Comparison of MTL (FDLINE, ULM) and RMTL (FDLM) line mod-els for the unbalanced fault test of Fig. 3.6.is slightly closer to the reference solution than FDLINE, but the errors for ULMincrease in the subsequent cycles. This is shown more clearly in Fig. 4.10. Thisis probably related to the less accurate determination by ULM of the propagationdelays.The errors for FDLM remain relatively low over the entire simulation time win-dow. It can be seen in Fig. 4.8 that there is a slight phase shift between the MTL andthe RMTL solutions that increases as time advances. This shift might be related tothe waves collocation issues discussed in this chapter. All simulations simulationsconverge when steady state is reached.Table 4.4 shows the maximum instantaneous errors for all voltages and cur-rents. The maximum relative errors for FDLINE and FDLM remained relatively950 2 4 6 8Time (ms)-500-400-300-200-1000v3 (kV)8 10 12 14 16Time (ms)-1500150300450v3 (kV)16 18 20 22 24Time (ms)-450-300-1500150v3 (kV)24 26 28 30 32Time (ms)-400-2000200400v3 (kV)DTFS: RMTLFDLMDTFS: MTLFDLINEULMFigure 4.9: Details of the first four peaks of v3 in Fig. 4.8a.even, while the maximum errors were larger for ULM for v1 and v4. The results forFDLINE and ULM were consistent with the previous results of [116].Table 4.4: Maximum errors for voltages and currents for FDLINE, ULM, andFDLM with respect to their reference DTFS solution.f FDLINE ULM FDLMv1 4.26% 10.50% 5.07%v3 7.18% 8.56% 6.67%v4 5.53% 10.48% 5.03%v5 8.19% 12.15% 5.42%i2 2.30% 2.33% 2.07%i6 3.10% 2.13% 2.04%960 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4Time (s)0246810Error of v3 (%)FDLINE ULM FDLMFigure 4.10: Errors for the voltage v3 for FDLINE and ULM versus DTFS forthe MTL equations and for FDLM versus DTFS for the RMTL equations.Computational efficiency of FDLINE, ULM, and FDLM are compared in Table4.5 in the simulation of unbalanced faults test of Fig. 3.6. For this comparison,operations count per time step is calculated for FDLINE and FDLM using (2.31) andfor ULM using (2.35). The fitting information were extracted from Microtran forFDLINE, PSCAD for ULM, and BAF algorithm of [1] for FDLM.Table 4.5: Comparison of operations count for FDLINE, ULM, and FDLM linemodels in the simulation of unbalanced faults test of Fig. 3.6.Models Zc or Yc e−γ` totalFDLINE 6×6+142 6×6+182 396FDLM 6×6+87 6×6+231 390ULM 6×6×11 6×6×41 1872As indicated in Table 4.5, ULM is about 4.7 times more expensive than FDLINEand FDLM in terms of computational cost. Even though, ULM requires fewer polesfor fitting compared to FDLINE and FDLM, the main computational burden is due tofitting the frequency dependence of the transformation matrices. Since FDLINE andFDLM use the same fitting routine, their computational cost are relatively similar.97Figure 4.11 show the results for the open circuit voltage v1 for the test of Fig.3.8 obtained with FDLM with T calculated at 1 Hz, 100 Hz, and 1 kHz comparedto the DTFS with T calculated at each frequency.As opposed to the results of Fig. 3.30 where the accuracy of FDLINE dependedon the chosen frequency, FDLM is a frequency-independent model.0 4 8 12 16 20Time (ms)-0.250.250.751.251.752.25v1 (kV)0.1 0.11 0.12 0.13 0.14 0.15Time (s)0.850.951.051.15v1 (kV)DTFS: T , Skin in C FDLM @ 1 Hz FDLM @ 60 Hz FDLM @ 1 kHzFigure 4.11: v1 in the step response test of Fig. 3.8. Comparison of RMTLmodels: DTFS versus FDLM (for FDLM T calculated at different fre-quencies).4.10 General ObservationsThe fundamental constraint for the MTL equations proposed in this chapter was thatconceptually the voltage and current waves must be collocated and travel togetherwith the same propagation function. Another factor for maintaining the physicalconsistency in the line equations is to make the ground resistance contributions toall elements of the Z matrix equal at each frequency. Based on this condition, theRMTL equations are proposed.98The diagonalization method to find real transformation matrices for multicon-ductor transmission lines described by the RMTL equations maintains the geometri-cal identity of the propagation modes by separating the inductance and capacitance(responsible for propagation) from the losses. The problem is then reduced to di-agonalizing a real matrix Lloop, which can be achieved with a simple eigenvectorsroutine. There are no mode switchings in this process. This process results ina real transformation matrix that is practically constant over the entire frequencyrange. Error analysis demonstrated the validity of this assumption for diagonaliza-tion of the individual Z and Y and the product ZY and YZ for an asymmetricalline configuration.The fact that it is possible to find a real constant transformation matrix to ac-curately diagonalize the line equations validates the physical assumption that thesematrices should be diagonalized using a real transformation matrix. It also explainswhy the FDLINE model of [69] which uses real constant TV and TI as an approxi-mation under the MTL equations, gives very accurate results, even for asymmetricalline configurations [116].A new FDLM transmission line model was proposed in this chapter to imple-ment the RMTL equations. For the proposed FDLM under the RMTL equations,using a real constant T for both voltages and currents is not an approximation butconceptually exact, within the limits of validity of these equations.Simulation results were used to compare the EMTP frequency dependent linemodels FDLINE and ULM, and the proposed FDLM model, against reference fre-quency domain solutions using the MTL and the RMTL equations. The simulationresults were very close for all models, with smaller overall errors for FDLM, fol-lowed by FDLINE. Overall, the results obtained with the proposed RMTL model arevery close to the solutions given by the classical MTL models. These results areexact within the first cycle of the transients which the most important cycle for thesetting of relays and protection devices.More fundamental analytical work is needed to incorporate additional symme-try considerations into the line equations. One of these considerations is the princi-ple of collocation of charge and current, which would require the re-collocation ofthe charge in the conductors and the ground when the current position is affectedby the skin effect. In the next section, we investigate this concept with the available99line equations.4.11 Further Investigations on the Wave CollocationConditionThe waves collocation condition of (4.12) leads naturally to the concept that theexternal inductance for the wave propagation Lloop and the external capacitancefor the wave propagation Cext should be symmetrical. However, while Lloop issubjected to skin effect, Cext is held constant at its electrostatic value by the fun-damental assumption of TEM propagation. This inconsistency, which violates theprinciple of collocation of charge and current, has also an influence on the colloca-tion of the voltage and current waves. In this chapter, we enforced the collocationof the voltage and current waves by requiring that ZY = YZ (wave collocation)but more research is needed to understand how the propagation equations can bemade consistent with the additional physical requirement of collocation of chargeand current.Skin effect in the capacitance has been described in the literature by as earlyas Wise [138], and later by Nakagawa [81] and other authors, however, it has notbeen incorporated in the EMTP line models.In this section, we explore the effect of considering the “skin effect” in thecapacitances on the RMTL equations proposed in this research knowing that wemay contradict the TEM condition.For consistency between the inductance and capacitance matrices, we use theLord Kelvin’s method of imaging described in Fig. 2.2 to calculate the self andmutual elements of the complex loop Maxwell’s potential coefficient matrix P′loop,P′loop11 =12piεln2(h1+ p)r1eq, P′loop12 =12piεlnD12′d12(4.33)Similarly to L′loop, the elements of C′loop comprise two parts: real and imagi-nary. The imaginary part in (4.33) will correspond to the transversal losses due tothe ground resistivity. The procedure to calculate the transversal losses are similarto the expressions (4.14). In Fig. 4.12, elements of earth transversal losses for theline of Fig. 3.5 are shown for a frequency range from 10−2 Hz to 106 Hz.10010-2 10-1 100 101 102 103 104 105 106Frequency (Hz)-10-2-10-4-10-6-10-8-10-10-10-12-10-14Earth transversal losses (S/km)Figure 4.12: Earth transversal losses given by the complex penetration depthformula for the double circuit vertical line of Fig. 3.5As observed in Fig. 4.12, earth transversal losses are very small quantities withnegative values, therefore, to maintain the passivity in the system these elementsare eliminated from the formulation.The matrix of capacitances Cloop is obtained by taking the inverse of the realpart of P′loop,Cloop = P−1loop (4.34)Figure 4.13 shows the elements of Lloop (valid for MTL and RMTL equations)in comparison with the elements of Cloop (investigated in this section).As observed in Fig. 4.13, in a frequency range of 10−2 Hz to 106 Hz, Lloopdecreases with frequency while Cloop slightly increases with frequency. The maxi-mum rate of change for Lloop is about 50% for diagonals and 85% for off-diagonals,whereas for Cloop is about 8% for diagonals and 50% for off-diagonals. With areasonable compromise that in Cloop the diagonals can be considered constant (tocomply with the TEM assumption) and the off-diagonals are less dominant withrespect to the diagonals when compared to the off-diagonals of Lloop, we use the101 10-210-110010110210310410510601234Lloop(mH/km)10-210-1100101102103104105106Frequency (Hz)-30369Cloop(nF/km)Diagonal Off-Diagonal Diagonal Off-Diagonal Figure 4.13: Lloop and Cloop for the line of Fig. 3.5.skin effect in the capacitances to further investigate the RMTL equations.First, we evaluate the effect of considering the skin effect in the capacitanceson the wave collocation condition (4.12). Figure 4.14 compares ZY(YZ)−1 for theline of Fig. 3.5 for two conditions: 1) no skin effect in the capacitances (Cext), 2)with skin effect in the capacitances (Cloop). The curves in Fig. 4.14 are obtainedusing the RMTL equations for both conditions.These results signify that considering the “skin effect” has a significant impacton the wave collocation condition by reducing the error of ZY(YZ)−1 to the max-imum of about 0.1% for the real part (barely visible in the plot) and 0.5% for theimaginary part.Next, we compare the effect of considering the skin effect in the capacitanceson the product of the inductance and capacitance matrices. This product is nor-malized by υ2c , where υc is the speed of light (3× 108 m/s). The diagonal andoff-diagonal elements of the product of inductance and capacitance matrices areplotted in Fig. 4.15.As observed in Fig. 4.15, the diagonal terms of the normalized LloopCloop are“One” and the off-diagonal terms are “Zero” for all frequencies. These results10210-2 10-1 100 101 102 103 104 105 1060.980.9911.011.02Real of ZY(YZ)-1 Diagonal No Skin in CSkin in C10-2 10-1 100 101 102 103 104 105 106-0.06-0.0300.030.06Real of ZY(YZ)-1 Off-Diagonal10-2 10-1 100 101 102 103 104 105 106Frequency (Hz)-0.03-0.02-0.0100.010.020.03Imaginary of ZY(YZ)-1All ElementsFigure 4.14: Validity of the collocation condition. If ZY = YZ, the prod-uct ZY(YZ)−1 should be one for the diagonals and zero for the off-diagonals.seem to comply with physics that the propagation function only depends on themedium (LC = µ0ε0 = 1/υ2). On the other hand, the elements of the normalizedLloopCext vary with frequency until very high frequencies that both inductance andcapacitance matrices are at their electrostatic position values.With this premise, we use the skin effect in the capacitances to further explorethe time domain transient simulation. For this analysis, we run the simulationcondition of of Fig. 4.8a for the case that skin effect is considered in capacitances.Figure 4.16 is used to test two conceptual in these experiments:i) MTL models (using two complex frequency dependent TV (ω) and TI(ω))10310-2 10-1 100 101 102 103 104 105 10611.52vc2 LloopCDiagonal Cloop Cext10-2 10-1 100 101 102 103 104 105 106Frequency (Hz)00.51vc2 LloopCOff-DiagonalFigure 4.15: The effect of skin effect in capacitances on the product LC forthe double circuit vertical line of Fig. 3.5.versus RMTL models (using a constant real T).ii) No skin effect in the capacitances (Fig. 4.16a) versus skin effect in thecapacitances (Fig 4.16b).Results of Fig. 4.16b suggest that including skin effect in the capacitances willfurther enforce the collocation wave condition as the MTL and RMTL models co-incide. As opposed to Fig. 4.16a, there is no shift between the MTL and RMTLmodels when skin effect is considered in capacitances. In other words, colloca-tion of charges and currents naturally contribute in the collocation of voltage andcurrent waves (with no enforcement).Figure 4.17 uses DTFS to compare the difference between the MTL equationswith no skin effect in the C versus RMTL equations with skin effect in C.1040 5 10 15 20 25 30 35 40 45 50-600-3000300600v3 (kV)(a)DTFS: RMTLFDLMDTFS: MTLFDLINEULM0 5 10 15 20 25 30 35 40 45 50Time (ms)-600-3000300600v3 (kV)(b)DTFS: RMTLFDLMDTFS: MTLFigure 4.16: Comparison of MTL and RMTL models for two conditions: a) noskin in C, b) with skin in C. v3 in the test of Fig. 3.6.0 5 10 15 20 25 30 35 40 45 50Time (ms)-600-3000300600v3 (kV)DTFS: MTL , No Skin in C DTFS: RMTL , Skin in CFigure 4.17: Comparison of MTL with Cext versus RMTL with Cloop. v3 in thetest of Fig. 3.6.105Even though the results presented in this section are very promising, the va-lidity of this premise still remains under the speculation since the equations thatwe used for this investigation are based on the TEM assumption that do not allowfrequency dependent capacitances.As observed in Fig. 4.17, considering skin effect in the capacitances has astronger influence in the simulation results. Further research is required to discoverthe equations that allow us to consider skin effect in the capacitances. Ultimatevalidation can be performed by comparing the results with a field test.106Chapter 5Summary of Contributions andFuture Works5.1 Conclusions and ContributionsThe main objective in this research was to develop a frequency dependent EMTPtransmission line model which, in addition to preserving high accuracy, numericalstability, high computational efficiency, and simplicity in implementation, coin-cides well with the fundamentals of electromagnetics. The second objective wasto verify the accuracy of the new line model with an accurate frequency domainreference solution. The steps taken to reach these objectives are as follows:Contribution 1The first step was to assess the accuracy of FDLINE model under asymmetricalline configurations. In Chapter 2, six different line configurations were simulatedunder asymmetrical short-circuit conditions using FDLINE and ULM. Three casesof three-phase lines were considered: single-circuit lines, double-circuit lines inthe same tower, and double-circuit lines in separate towers. Open-circuit voltagesand short-circuit currents were compared. The accuracy of time domain simula-tions were assessed using a conventional frequency-domain NLT solution. For allcases, the FDLINE model gave similar results to the ULM model and both mod-els gave good results when compared to the reference conventional NLT solution107for single-circuit lines. These results indicate that, contrary to traditional belief, aconstant transformation matrix model like FDLINE is capable of representing multi-circuit asymmetrical line configurations. With this validation, FDLINE was used inChapter 4 to develop a new line model based on the RMTL equations.In the simulation of double circuit lines in Chapter 2, we observed that theresults obtained with NLT solution did not match with the EMTP line models. Thishappened because the choices of damping factor and windowing function for NLTwere not suitable for those simulations.Contribution 2In Chapter 3, a simpler methodology based on the DTFS formulation was pre-sented which provides very accurate solutions without the need, unlike the NLT, toadjust additional parameters. The main premise for the application of DTFS is thatthe user specifies a time window width and a frequency window width as startingpoints for the solution. The time window width is determined by the system timeconstants and the frequency window width is determined by the desired maximumfrequency of the transient. Guidelines were provided for this set up. This approachguarantees a one-to-one matching between time points and frequency points andavoids the typical numerical errors of frequency domain solutions such as aliasingand Gibbs. The accuracy of the DTFS was tested for a variety of simulations in-cluding circuits with lumped elements and frequency dependent transmission lines.Contribution 3In Chapter 4, a new EMTP transmission line “FDLM” was developed that physi-cally satisfies the symmetry for the propagation of voltage waves and current wavesreferred to as “collocation wave condition”. This condition was reached by usinga single transformation matrix as a constraint to diagonalize the wave functions.The validity of this approximation was tested mathematically with high accuracy.The nature of real and constant of the transformation matrix within wide rangeof switching transients frequency further facilitated the diagonalization. Based onthese considerations, the RMTL equations were proposed which were further usedin the development of the FDLM model. Unifying the slightly different elements ofthe ground return resistance matrix given by the method of image or Carson’s for-108mula resulted in physical consistency and higher symmetry in the RMTL equations.Simulation results validated the accuracy of the FDLM compared to FDLINE,ULM, and newly developed DTFS reference solution for the double circuit verticalline. In terms of computational cost, FDLM is as efficient as FDLINE while usinga constant real transformation matrix is not an approximation but is exact over theentire frequency range for switching transients.5.2 Future WorksThe accuracy and functionality of the proposed FDLM model can be more justifiedwith performing the following studies:Study 1: Inclusion of ground wires and bundled conductors in the modelGround wires are steel conductors used to protect the overhead transmissionlines from the lightning. These conductors can be grounded at both ends of towersand they are eliminated from the matrix of impedance and admittance by using asfor example a Kron reduction when calculating line parameters. Since the radiusand DC resistance of ground wires are different from the conductors, the effect ofmatrix reduction may cause some asymmetries in the model. This is also the casefor bundled conductors when matrix reduction is used to propagate the effect ofbundling on other phases.When ground wires are segmented, they may contribute to the admittance ma-trix creating another source of asymmetry between the impedance and admittancematrices. These aspects should be investigated in the future work.Study 2: Test the RMTL model for underground cablesAs addressed in [71], in the MTL model the transformation matrices will be-come strongly frequency dependent in the case of the underground cables due tothe short distances from the conductors to the ground. As further extension to thisresearch, it would be interesting to reformulate the cable equations based on theproposed RMTL model and investigate the behaviour of the transformation matri-ces and the validity of the collocation wave condition.109Study 3: Test the RMTL model for a field testIn this thesis, the proposed FDLM passed the preliminary verifications throughcomparisons with the DTFS reference frequency domain solution for the same sys-tem topology and parameters. Since this model belongs to the new class of lineequations where voltage and current waves are collocated, the final verification willbe achieved by comparing the transient solution of the new line model with a fieldtest with an untransposed asymmetrical line configuration. This requires to accessto accurate measurement of the field test voltages and currents, line configuration,conductor type, equivalent source representation, and earth resistance.Study 4: Collocation of charges and currentsOne consideration to incorporate additional symmetry into the line equationsis the principle of collocation of charge and current, which would require the re-collocation of the charge in the conductors and the ground when the current po-sition is affected by the skin effect. In this research, we investigated this conceptwithin certain validity of TEM assumption. More fundamental research is requiredto find the line equations that can allow the use of skin effect in capacitances.110Bibliography[1] Improvements on the JMARTI transmission Line Model. DCG-EPRICoordination Group, 1983. → pages 2, 9, 23, 93, 97[2] Microtran Power Systems Analysis Corporation, MicroTran-TransientAnalysis Program Reference Manual, Vancouver, 2014. → pages 32, 72[3] H. Abdollahzadeh, M. Jazaeri, and A. Tavighi. A new fast-convergedestimation approach for dynamic voltage restorer (dvr) to compensatevoltage sags in waveform distortion conditions. International Journal ofElectrical Power & Energy Systems, 54:598–609, 2014. → pages 49[4] F. L. Alvarado and R. Betancourt. An accurate closed-form approximationfor ground return impedance calculations. Proceedings of the IEEE, 71(2):279–280, February 1983. → pages 4[5] A. Ametani. The application of the fast fourier transform to electricaltransient phenomena. International Journal of Electrical EngineeringEducation, 10(4):277–287, 1973. → pages 7[6] A. Ametani. Stratified earth effects on wave propagation frequencydependent parameters. IEEE Transactions on Power Apparatus andSystems, PAS-93(5):1233–1239, September 1974. → pages 4[7] A. Ametani and K. Imanishi. Development of exponential fouriertransform and its application to electrical transients. Proceedings of theInstitution of Electrical Engineers, 126(1):51–56, January 1979. → pages 7[8] A. Ametani and R. Schinzinger. Equations for surge impedance andpropagation constant of transmission lines above stratified earth. IEEETransactions on Power Apparatus and Systems, 95(3):773–781, May 1976.→ pages 4111[9] A. Arismunandar. Capacitive correction factors for transmission lines toinclude finite conductivity and dielectric constant of the earth. IEEETranssactions on Power Apparatus and Systems, pages 436–456, 1963. →pages 5[10] L. J. B. Bergeron. Water hammer in hydraulics and wave surges inelectricity. Wiley, 1961. → pages 8[11] J.-F. Blais, M. Cimmino, A. Ross, and D. Granger. Suppression of timealiasing in the solution of the equations of motion of an impacted beamwith partial constrained layer damping. Journal of Sound and Vibration,326(3–5):870–882, October 2009. → pages 55[12] A. Budner. Introduction of frequency-dependent line parameters into anelectromagnetic transients program. IEEE Transactions on PowerApparatus and Systems, PAS-89(1):88–97, January 1970. → pages 8[13] J. R. Carson. Wave propagation in overhead wires with ground return. BellSystem Technical Journal, 5(4):539–554, October 1926. → pages 3, 82[14] F. Castellanos and J. R. Martı´. Phase-domain multi-phase transmission linemodels. International Conference on Power System Transients (IPST 95),Lisbon, Portugal, pages 17–22, 3–7 September 1995. → pages 11[15] F. Castellanos and J. R. Martı´. Full frequency-dependent phase-domaintransmission line model. IEEE Transactions on Power Systems, 12(3):1331–1339, August 1997. → pages 10[16] A. Chinea and S. Grivet-Talocia. A passivity enforcement scheme fordelay-based transmission line macromodels. IEEE Microwave andWireless Components Letters, 17(8):562–564, August 2007. → pages 11[17] A. I. Chrysochos, T. A. Papadopoulos, and G. K. Papagiannis. Enhancingthe frequency-domain calculation of transients in multiconductor powertransmission lines. Electric Power Systems Research, 122:56–64, 2015. →pages 7[18] S. Cristina and M. D’Amore. Dft-based procedure for transmission-linetransient computation. IEE Proceedings C-Generation, Transmission andDistribution, 134(2):138–144, March 1987. → pages 7[19] C. G. Cullen. Matrices and linear transformations. Courier Corporation,2012. → pages 10112[20] M. D’Amore and M. S. Sarto. A new formulation of lossy ground returnparameters for transient analysis of multiconductor dissipative lines. IEEETransactions on Power Delivery, 12(1):303–314, January 1997. → pages 6[21] M. D’Amore and M. S. Sarto. Simulation models of a dissipativetransmission line above a lossy ground for a wide-frequency range. i. singleconductor configuration. IEEE Transactions on ElectromagneticCompatibility, 38(2):127–138, May 1996. → pages 6[22] M. D’Amore and M. S. Sarto. Simulation models of a dissipativetransmission line above a lossy ground for a wide-frequency range. ii.multiconductor configuration. IEEE Transactions on ElectromagneticCompatibility, 38(2):139–149, May 1996. → pages 6[23] S. J. Day, N. Mullineux, and J. R. Reed. Developments in obtainingtransient response using fourier transforms part i: Gibbs phenomena andfourier integrals. International Journal of Electrical EngineeringEducation, 3(4):501–506, 1965. → pages 7[24] S. J. Day, N. Mullineux, and J. R. Reed. Developments in obtainingtransient response using fourier transforms part ii: Use of the modifiedfourier transform. International Journal of Electrical EngineeringEducation, 4(1):31–40, 1966. → pages 7[25] H. M. J. De Silva, A. M. Gole, J. E. Nordstrom, and L. M. Wedepohl.Robust passivity enforcement scheme for time-domain simulation ofmulti-conductor transmission lines and cables. IEEE Transactions onPower Delivery, 25(2):930–938, April 2010. → pages 11[26] A. Deri, G. Tevan, A. Semlyen, and A. Castanheira. The complex groundreturn plane a simplified model for homogeneous and multi-layer earthreturn. IEEE Transactions on Power Apparatus and Systems, PAS-100(8):3686–3693, August 1981. → pages 3, 5, 16[27] D. Deschrijver, B. Gustavsen, and T. Dhaene. Causality preservingpassivity enforcement for traveling-wave-type transmission-line models.IEEE Transactions on Power Delivery, 24(4):2461–2462, October 2009.→ pages 11[28] H. W. Dommel. Electromagnetic transients program (EMTP) theory book.Bonneville Power Administration, 1986. → pages 2, 41113[29] H. W. Dommel. Digital computer solution of electromagnetic transients insingle and multiphase networks. IEEE Transactions on Power Apparatusand Systems, PAS-88(4):388–399, April 1969. → pages 2[30] H. W. Dommel and W. S. Meyer. Computation of electromagnetictransients. Proceedings of the IEEE, 62(7):983–993, July 1974. → pages 2[31] C. Dubanton. Calcul approche´ des parametres primaires et secondairesdune ligne de transport. EDF Bulletin de la Direction des Etudes etRecherches, 1:53–62, 1969. → pages 3[32] C. Dufour, H. Le-Huy, J. C. Soumagne, and A. El Hakimi. Real-timesimulation of power transmission lines using martı´ model with optimalfitting on dual-dsp card. IEEE Transactions on Power Delivery, 11(1):412–419, January 1996. → pages 9[33] A. E. Efthymiadis and L. M. Wedepohl. Propagation characteristics ofinfinitely-long single-conductor lines by the complete field solutionmethod. Proceedings of the Institution of Electrical Engineers, 125(6):511–517, June 1978. → pages 6[34] E. P. R. I. (EPRI). Transmission-line reference book. 345 kv and above.Technical report, 2nd Edition, General Electric Co., Pittsfield, MA (USA).Large Transformer Div.; General Electric Co., Schenectady, NY (USA).Electric Utility Systems Engineering Dept., 1982. → pages 28, 51, 94[35] J. C. Escamilla, P. Moreno, and P. Go´mez. New model for overhead lossymulticonductor transmission lines. IET Generation, Transmission &Distribution, 7(11):1185–1193, 2013. → pages 10[36] J. A. B. Faria and J. F. B. Da Silva. Wave propagation in polyphasetransmission lines a general solution to include cases where ordinary modaltheory fails. IEEE Transactions on Power Delivery, 1(2):182–189, April1986. → pages 9[37] J. A. B. Faria and J. F. B. da Silva. The effect of randomly earthed groundwires on plc transmission-a simulation experiment. IEEE Transactions onPower Delivery, 5(4):1669–1677, November 1990. → pages 10[38] J. B. Faria. Overhead three-phase transmission lines-nondiagonalizablesituations. IEEE Transactions on Power Delivery, 3(4):1348–1355,October 1988. → pages 9114[39] A. B. Fernandes, W. L. Neves, E. G. Costa, and M. N. Cavalcanti.Transmission line shunt conductance from measurements. IEEETransactions on Power Delivery, 19(2):722–728, 2004. → pages 16, 51[40] R. H. Galloway, W. B. Shorrocks, and L. M. Wedepohl. Calculation ofelectrical parameters for short and long polyphase transmission lines.Proceedings of the Institution of Electrical Engineers, 111(12):2051–2059,1964. → pages 5[41] C. Gary. Approche comple`te de la propagation multifilaire en hautefre´quence par utilisation des matrices complexes. EDF Bulletin de laDirection des Etudes et Recherches, 3(4):5–20, 1976. → pages 3[42] P. Go´mez and J. C. Escamilla. Frequency domain modeling of nonuniformmulticonductor lines excited by indirect lightning. International Journal ofElectrical Power & Energy Systems, 45(1):420–426, 2013. → pages 7[43] P. Go´mez and F. A. Uribe. The numerical laplace transform: An accuratetechnique for analyzing electromagnetic transients on power systemdevices. International Journal of Electrical Power & Energy Systems, 31(2–3):116–123, February–March 2009. → pages 7, 26, 27[44] P. Go´mez, P. Moreno, and J. L. Naredo. Frequency-domain transientanalysis of nonuniform lines with incident field excitation. IEEETransactions on Power Delivery, 20(3):2273–2280, July 2005. → pages 7[45] B. Gustavsen. Validation of frequency-dependent transmission line models.IEEE Transactions on Power Delivery, 20(2):925–933, April 2005. →pages 7[46] B. Gustavsen. Improving the pole relocating properties of vector fitting.IEEE Transactions on Power Delivery, 21(3):1587–1592, July 2006. →pages 11, 25[47] B. Gustavsen. Avoiding numerical instabilities in the universal line modelby a two-segment interpolation scheme. IEEE Transactions on PowerDelivery, 28(3):1643–1651, July 2013. → pages 11[48] B. Gustavsen. Passivity enforcement for transmission line models based onthe method of characteristics. IEEE Transactions on Power Delivery, 23(4):2286–2293, October 2008. → pages 11115[49] B. Gustavsen. Modal domain-based modeling of parallel transmission lineswith emphasis on accurate representation of mutual coupling effects. IEEETransactions on Power Delivery, 27(4):2159–2167, October 2012. →pages 10, 28[50] B. Gustavsen and J. Nordstrom. Pole identification for the universal linemodel based on trace fitting. IEEE Transactions on Power Delivery, 23(1):472–479, January 2008. → pages 11[51] B. Gustavsen and A. Semlyen. Combined phase and modal domaincalculation of transmission line transients based on vector fitting. IEEETransactions on Power Delivery, 13(2):596–604, April 1998. → pages 11[52] B. Gustavsen and A. Semlyen. Simulation of transmission line transientsusing vector fitting and modal decomposition. IEEE Transactions onPower Delivery, 13(2):605–614, April 1998. → pages 11[53] B. Gustavsen and A. Semlyen. Enforcing passivity for admittance matricesapproximated by rational functions. IEEE Transactions on Power Systems,16(1):97–104, February 2001. → pages 11[54] B. Gustavsen and A. Semlyen. Calculation of transmission line transientsusing polar decomposition. IEEE Transactions on Power Delivery, 13(3):855–862, July 1998. → pages 11[55] B. Gustavsen and A. Semlyen. Rational approximation of frequencydomain responses by vector fitting. IEEE Transactions on Power Delivery,14(3):1052–1061, July 1999. → pages 2, 11, 23[56] D. E. Hedman. Propagation on overhead transmission lines: I-theory ofmodal analysis, ii-earth conduction effects and practical results. IEEETransactions on Power Apparatus and Systems, 84(3):200–211, March1965. → pages 5, 8, 20, 86[57] W. Hendrickx and T. Dhaene. A discussion of “rational approximation offrequency domain responses by vector fitting”. IEEE Transactions onPower Systems, 21(1):441–443, February 2006. → pages 24[58] L. Hofmann. Series expansions for line series impedances consideringdifferent specific resistances, magnetic permeabilities, and dielectricpermittivities of conductors, air, and ground. IEEE Transactions on PowerDelivery, 18(2):564–570, April 2003. → pages 3116[59] K. Iwamoto. Use of travelling waves on the measurement of earthresistivity. Journal of Institution of Electrical Engineers of Japan, 78:1038–1049, 1958. → pages 4[60] K. Iwamoto. Wide time-domain calculation of energizing surges inmulti-conductor transmission line by a method of logarithmic fouriertransformation. IEEE Transactions on Power Apparatus and Systems,PAS-100(4):1788–1795, April 1981. → pages 7[61] H. Keshtkar, S. K. Solanki, and J. M. Solanki. Improving the accuracy ofimpedance calculation for distribution power system. IEEE Transactionson Power Delivery, 29(2):570–579, April 2014. → pages 4[62] H. Kikuchi. Wave propagation along an infinite wire above ground at highfrequencies. Electrotechnical Journal of Japan, 2:73–78, 1956. → pages 5[63] I. Kocar, J. Mahseredjian, and G. Olivier. Improvement of numericalstability for the computation of transients in lines and cables. IEEETransactions on Power Delivery, 25(2):1104–1111, April 2010. → pages11[64] O. R. Lean˜os, J. L. Naredo, and P. Moreno. Assessment of approximateformulas for calculating overhead-line earth-impedances. 40th NorthAmerican Power Symposium, pages 1–6, 28–30 September 2008. → pages4[65] F. J. Marcano. Modelling of transmission lines using idempotentdecomposition. PhD thesis, University of British Columbia, August 1996.→ pages 11[66] F. J. Marcano and J. R. Martı´. Idempotent line model: Case studies.International Conference on Power System Transients (IPST 97), Seattle,Washington, pages 67–72, 22–26 June 1997. → pages 11, 28[67] J. R. Martı´. Digital signal processing from a simulation point of view. classnotes for elec 466, 1999. → pages 41[68] J. R. Martı´. The problem of frequency dependence in transmission linemodelling. PhD thesis, University of British Columbia, April 1981. →pages 2, 8, 9, 23[69] J. R. Martı´. Accurate modelling of frequency-dependent transmission linesin electromagnetic transient simulations. IEEE Transactions on Power117Apparatus and Systems, PAS-101(1):147–157, January 1982. → pages 2,9, 22, 23, 92, 99[70] J. R. Martı´, H. W. Dommel, L. Martı´, and V. Brandwajn. Approximatetransformation matrices for unbalanced transmission lines. 9th PowerSystem Computation Conference (PSCC 87), Butterworths, London, pages416–422, August 1987. → pages 9[71] L. Martı´. Simulation of transients in underground cables withfrequency-dependent modal transformation matrices. IEEE Transactionson Power Delivery, 3(3):1099–1110, July 1988. → pages 2, 10, 109[72] L. Martı´, R. Y. Brierley, and T. E. Grainger. Analysis of electromagnetictransients in cross-bounded cable systems using frequency dependent cablemodels. International Conference on Power System Transients (IPST 95),Lisbon, Portugal, pages 5–10, 3–7 September 1995. → pages 10[73] W. S. Meyer and H. W. Dommel. Numerical modelling offrequency-dependent transmission-line parameters in an electromagnetictransients program. IEEE Transactions on Power Apparatus and Systems,PAS-93(5):1401–1409, September 1974. → pages 8[74] A. Morched, B. Gustavsen, and M. Tartibi. A universal model for accuratecalculation of electromagnetic transients on overhead lines andunderground cables. IEEE Transactions on Power Delivery, 14(3):1032–1038, July 1999. → pages 2, 11, 24[75] P. Moreno and A. Ramirez. Implementation of the numerical laplacetransform: a review task force on frequency domain methods for emtstudies, working group on modeling and analysis of system transients usingdigital simulation, general systems subcommittee, ieee power engineeringsociety. IEEE Transactions on Power Delivery, 23(4):2599–2609, October2008. → pages 7, 25, 26, 55[76] P. Moreno, P. Go´mez, J. L. Naredo, and J. L. Guardado. Frequency domaintransient analysis of electrical networks including non-linear conditions.International Journal of Electrical Power & Energy Systems, 27(2):139–146, February 2005. → pages 7[77] P. Moreno, R. de la Rosa, and J. L. Naredo. Frequency domaincomputation of transmission line closing transients. IEEE Transactions onPower Delivery, 6(1):275–281, January 1991. → pages 7118[78] N. Mullineux, J. R. Reed, and L. M. Wedepohl. Calculation of electricalparameters for short and long polyphase transmission lines. Proceedings ofthe Institution of Electrical Engineers, 112(4):741–742, April 1965. →pages 4[79] N. Nagaoka and A. Ametani. A development of a generalizedfrequency-domain transient program-ftp. IEEE Transactions on PowerDelivery, 3(4):1996–2004, October 1988. → pages 7[80] M. Nakagawa. Further studies on wave propagation along overheadtransmission lines: Effects of admittance correction. IEEE Transactions onPower Apparatus and Systems, PAS-100(7):3626–3633, July 1981. →pages 5, 6[81] M. Nakagawa. Admittance correction effects of a single overhead line.IEEE Transactions on Power Apparatus and Systems, PAS-100(3):1154–1161, March 1981. → pages 5, 6, 100[82] M. Nakagawa and K. Iwamoto. Earth-return impedance for the multi-layercase. IEEE Transactions on Power Apparatus and Systems, 95(2):671–676,April 1976. → pages 4[83] M. Nakagawa, A. Ametani, and K. Iwamoto. Further studies on wavepropagation in overhead lines with earth return: impedance of stratifiedearth. Proceedings of the Institution of Electrical Engineers, 120(12):1521–1528, December 1973. → pages 4[84] J. L. Naredo, A. C. Soudack, and J. R. Martı´. Simulation of transients ontransmission lines with corona via the method of characteristics. IEEProceedings-Generation, Transmission and Distribution, 142(1):81–87,January 1995. → pages 10[85] H. V. Nguyen, H. W. Dommel, and J. R. Martı´. Direct phase-domainmodelling of frequency-dependent overhead transmission lines. IEEETransactions on Power Delivery, 12(3):1335–1342, July 1997. → pages 10[86] T. Noda. Numerical techniques for accurate evaluation of overhead lineand underground cable constants. IEEJ Transactions on Electrical andElectronic Engineering, 3(5):549–559, August 2008. → pages 4[87] T. Noda. A double logarithmic approximation of carson’s ground-returnimpedance. IEEE Transactions on Power Delivery, 21(1):472–479, January2005. → pages 4119[88] T. Noda, N. Nagaoka, and A. Ametani. Phase domain modeling offrequency-dependent transmission lines by means of an arma model. IEEETransactions on Power Delivery, 11(1):401–411, January 1996. → pages10, 25[89] R. Nuricumbo-Guillen, P. Gomez, F. P. Espino-Cortes, and F. A. Uribe.Accurate computation of transient profiles along multiconductortransmission systems by means of the numerical laplace transform. IEEETransactions on Power Delivery, 29(5):2385–2393, 2014. → pages 7[90] T. A. Papadopoulos, G. K. Papagiannis, and D. A. Labridis. Wavepropagation characteristics of overhead conductors above imperfectstratified earth for a wide frequency range. IEEE Transactions onMagnetics, 45(3):1064–1067, March 2009. → pages 6[91] G. K. Papagiannis, D. A. Tsiamitros, D. P. Labridis, and P. S. Dokopoulos.A systematic approach to the evaluation of the influence of multilayeredearth on overhead power transmission lines. IEEE Transactions on PowerDelivery, 20(4):2594–2601, October 2005. → pages 4[92] A. Papoulis. The Fourier Integral and its Applications. McGraw-Hill, NewYork, 1962. → pages 23[93] C. R. Paul. Analysis of multiconductor transmission lines. 2nd Edition,John Wiley & Sons, 2008. → pages 19, 79, 80[94] P. Pettersson. Image representation of wave propagation on wires above, onand under ground. IEEE Transactions on Power Delivery, 9(2):1049–1055,April 1994. → pages 6[95] P. Pettersson. Propagation of waves on a wire above a lossyground-different formulations with approximations. IEEE Transactions onPower Delivery, 14(3):1173–1180, July 1999. → pages 6[96] M. Pizarro and R. Eriksson. Modeling of the ground mode of transmissionlines in time domain simulations. 7th international symposium on highvoltage engineering, (ISH 91), pages 179–182, 26–30 August 1991. →pages 4[97] F. Pollaczek. On the induction effects of a single phase ac line. ElectrischeNachricten Technik, 4:18–30, April 1927. → pages 3120[98] F. Pollaczek. On the field produced by an infinitely long wire carryingalternating current. Electrische Nachricten Technik, 3:339–359, March1926. → pages 3[99] J. G. Proakis and D. G. Manolakis. Digital Signal Processing: Principles,Algotithms, and Applications. 3rd Edition, Prentice Hall, 1996. → pages 49[100] A. Ramirez and F. Uribe. A broad range algorithm for the evaluation ofcarson’s integral. IEEE Transactions on Power Delivery, 22(2):1188–1193,April 2007. → pages 4[101] A. Ramirez, P. Gomez, P. Moreno, and A. Gutierrez. Frequency domainanalysis of electromagnetic transients through the numerical laplacetransforms. In Power Engineering Society General Meeting (PES 2004),pages 1136–1139. IEEE, 2004. → pages 7[102] A. Ramı´rez, J. L. Naredo, and P. Moreno. Full frequency-dependent linemodel for electromagnetic transient simulation including lumped anddistributed sources. IEEE Transactions on Power Delivery, 20(1):292–299,January 2005. → pages 10[103] C. Sanathanan and J. Koerner. Transfer function synthesis as a ratio of twocomplex polynomials. IEEE Transactions on Automatic Control, 8(1):56–58, January 1963. → pages 24[104] A. Semlyen. Some frequency domain aspects of wave propagation onnonuniform lines. IEEE Transactions on Power Delivery, 18(1):315–322,2003. → pages 7[105] A. Semlyen. Correspondence: Approximation to carson’s loss formulae.Canadian Electrical Engineering Journal, 6(2):30–31, April 1981. →pages 3[106] A. Semlyen. Ground return parameters of transmission lines an asymptoticanalysis for very high frequencies. IEEE Transactions on Power Apparatusand Systems, PAS-100(3):1031–1038, March 1981. → pages 3[107] A. Semlyen and A. Dabuleanu. Fast and accurate switching transientcalculations on transmission lines with ground return using recursiveconvolutions. IEEE Transactions on Power Apparatus and Systems, 94(2):561–571, March 1975. → pages 9121[108] A. Semlyen and B. Gustavsen. Phase-domain transmission-line modelingwith enforcement of symmetry via the propagated characteristic admittancematrix. IEEE Transactions on Power Delivery, 27(2):626–631, April 2012.→ pages 11, 80[109] A. Semlyen and A. Ramirez. Direct frequency domain computation oftransmission line transients due to switching operations. IEEE Transactionson Power Delivery, 23(4):2255–2261, October 2008. → pages 7[110] A. Semlyen and D. Shirmohammadi. Calculation of induction andmagnetic field effects of three phase overhead lines above homogeneousearth. IEEE Transactions on Power Apparatus and Systems, PAS-101(8):2747–2754, September 1982. → pages 3[111] J. K. Snelson. Propagation of travelling waves on transmissionlines-frequency dependent parameters. IEEE Transactions on PowerApparatus and Systems, PAS-91(1):85–91, January 1972. → pages 8[112] A. Sommerfeld. Propagation of waves in wireless telegraphy. Ann. derPhys, 28(3):665–736, March 1909. → pages 3[113] E. D. Sunde. Earth conduction effects in transmission systems. DoverPublications Inc., 1949. → pages 4, 5[114] M. C. Tavares, J. Pissolato, and C. M. Portela. Mode domain multiphasetransmission line model-use in transient studies. IEEE Transactions onPower Delivery, 14(4):1533–1544, 1999. → pages 10[115] A. Tavighi, S. Soleymani, J. Martı´, and H. Abdollahzadeh. Swift detectionof power transformer inrush current based on the windowed-adaptive linearcombiner estimation algorithm. → pages 49[116] A. Tavighi, J. R. Martı´, and J. A. Gutierrez-Robles. Comparison of thefdline and ulm frequency dependent emtp line models with a referencelaplace solution. International Conference on Power System Transients(IPST 15), Cavtat, Croatia, 15–18 June 2015. → pages 50, 51, 94, 96, 99[117] A. Tavighi, H. Abdollahzadeh, and J. Martı´. Fast response dvr controlstrategy design to compensate unbalanced voltage sags and swells indistribution systems. In Power and Energy Society General Meeting (PES),2013 IEEE, pages 1–5. IEEE, 2013. → pages 49122[118] T. Theodoulidis. On the closed-form expression of carson’s integral.Periodica Polytechnica. Electrical Engineering and Computer Science, 59(1):26–29, 2015. → pages 82[119] D. A. Tsiamitros, G. K. Papagiannis, and P. S. Dokopoulos. Homogenousearth approximation of two-layer earth structures: An equivalent resistivityapproach. IEEE Transactions on Power Delivery, 22(1):658–666, January2007. → pages 4[120] D. A. Tsiamitros, G. K. Papagiannis, and P. S. Dokopoulos. Earth returnimpedances of conductor arrangements in multilayer soils–part ii:numerical results. IEEE Transactions on Power Delivery, 23(4):2401–2408, November 2008. → pages 4[121] D. A. Tsiamitros, G. K. Papagiannis, and P. S. Dokopoulos. Earth returnimpedances of conductor arrangements in multilayer soils–part i:theoretical model. IEEE Transactions on Power Delivery, 23(4):2392–2400, October 2008. → pages 4[122] F. A. Uribe, J. L. Naredo, P. Moreno, and L. Guardado. Electromagnetictransients in underground transmission systems through the numericallaplace transform. International Journal of Electrical Power & EnergySystems, 24(3):215–221, 2002. → pages 7[123] S. Vujevic´ and D. Lovric´. Inverse continuous numerical fourier transformfor transient analysis of electromagnetic phenomena. IEEE Transactionson Electromagnetic Compatibility, 57(5):1149–1154, October 2015. →pages 41[124] J. R. Wait. Electromagnetic wave theory. Harper & Row, 1985. → pages90[125] J. R. Wait and K. P. Spies. On the image representation of the quasi-staticfields of a line current source above the ground. Canadian Journal ofPhysics, 47(23):2731–2733, May 1969. → pages 6[126] L. M. Wedepohl. Frequency domain analysis of wave propagation inmulticonductor transmission systems. Lecture Notes, University of BritishColumbia, Department of Electrical and Computer Engineering, Canada,1993. → pages 10, 19, 79[127] L. M. Wedepohl. Power system transients: Errors incurred in the numericalinversion of the laplace transform. 26th Midwest Symposium on Circuitsand Systems, (MSCS 83), pages 174–178, August 1983. → pages 7, 26123[128] L. M. Wedepohl. Application of matrix methods to the solution oftravelling-wave phenomena in polyphase systems. Proceedings of theInstitution of Electrical Engineers, 110(12):2200–2212, December 1963.→ pages 8, 20, 86[129] L. M. Wedepohl and A. E. Efthymiadis. Wave propagation in transmissionlines over lossy ground: a new, complete field solution. Proceedings of theInstitution of Electrical Engineers, 125(6):505–510, June 1978. → pages 6[130] L. M. Wedepohl and C. S. Indulkar. Switching overvoltages in shortcrossbonded cable systems using the fourier transform. Proceedings of theInstitution of Electrical Engineers, 122(11):1217–1221, November 1975.→ pages 7[131] L. M. Wedepohl and S. E. T. Mohamed. Transient analysis ofmulticonductor transmission lines with special reference to nonlinearproblems. Proceedings of the Institution of Electrical Engineers, 117(5):979–988, May 1970. → pages 7[132] L. M. Wedepohl and S. E. T. Mohamed. Multiconductor transmission lines.theory of natural modes and fourier integral applied to transient analysis.Proceedings of the Institution of Electrical Engineers, 116(9):1553–1563,September 1969. → pages 7[133] L. M. Wedepohl and R. G. Wasley. Wave propagation in multiconductoroverhead lines. calculation of series impedance for multilayer earth.Proceedings of the Institution of Electrical Engineers, 113(4):627–632,April 1966. → pages 5[134] L. M. Wedepohl and D. J. Wilcox. Transient analysis of undergroundpower-transmission systems. system-model and wave-propagationcharacteristics. Proceedings of the Institution of Electrical Engineers, 120(2):253–260, February 1973. → pages 7[135] L. M. Wedepohl, H. V. Nguyen, and G. D. Irwin. Frequency-dependenttransformation matrices for untransposed transmission lines usingnewton-raphson method. IEEE Transactions on Power Systems, 11(3):1538–1546, August 1996. → pages 21, 68, 87[136] D. J. Wilcox. Numerical laplace transformation and inversion.International Journal of Electrical Engineering Education, 15(3):247–265,July 1978. → pages 7, 26124[137] W. H. Wise. Propagation of high-frequency currents in ground returncircuits. Proceedings of the Institute of Radio Engineers, 22(4):522–527,April 1934. → pages 5[138] W. H. Wise. Potential coefficients for ground return circuits. Bell SystemTechnical Journal, 27(2):365–371, April 1948. → pages 5, 100[139] W. H. Wise. Effect of ground permeability on ground return circuits. BellSystem Technical Journal, 10(3):472–484, July 1931. → pages 5[140] L. Wu, P. A. A. F. Wouters, and F. E. Steennis. Frequency-domain transientanalysis in double-circuit mixed hv overhead line–cable connectionincluding cross-bonding. International Transactions on Electrical EnergySystems, 26(7):1408–1426, October 2015. → pages 41[141] T.-C. Yu and J. R. Martı´. A robust phase-coordinates frequency-dependentunderground cable model (zcable) for the emtp. IEEE Transactions onPower Delivery, 18(1):189–194, January 2003. → pages 10125
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A frequency-dependent multiconductor transmission line...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A frequency-dependent multiconductor transmission line model with collocated voltage and current propagation Tavighi, Arash 2017
pdf
Page Metadata
Item Metadata
Title | A frequency-dependent multiconductor transmission line model with collocated voltage and current propagation |
Creator |
Tavighi, Arash |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | This research contributes to developing a time domain and a frequency domain formulations to solve electromagnetic transients in power system with multiconductor overhead transmission lines. The time domain solution introduces a frequency dependent transmission line model “FDLM”. For the development of the FDLM a fundamental constraint is added to the classical line equations to maintain the symmetry between electric and magnetic fields. As a result, voltage waves and current waves travel together and the characteristic impedance remains uniform along the line. With this premise, a constant real transformation matrix can be obtained to diagonalize the line functions with high accuracy. This feature can greatly facilitate the line modelling as opposed to the existing line models which require complex frequency dependent transformation matrices for their diagonalization. The use of a single constant real transformation matrix for the voltage and current waves which is exact over the frequency range enables FDLM to provide higher accuracy and numerical efficiency than the existing line models while it complies with the physical system. The accuracy of the FDLM is assessed through comparisons with a newly developed Discrete Time Fourier Series frequency domain solution. This methodology is based on the correct specification of the time window and frequency window widths. Guidelines are provided for this set up which avoids the typical Gibbs and aliasing errors related to the classical frequency domain solutions. The proposed frequency domain solution is simpler to implement than the most commonly used numerical Laplace transform solution while it does not require further considerations to use damping factors or windowing functions. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-03-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0343065 |
URI | http://hdl.handle.net/2429/60791 |
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 |
GraduationDate | 2017-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_may_tavighi_arash.pdf [ 4.6MB ]
- Metadata
- JSON: 24-1.0343065.json
- JSON-LD: 24-1.0343065-ld.json
- RDF/XML (Pretty): 24-1.0343065-rdf.xml
- RDF/JSON: 24-1.0343065-rdf.json
- Turtle: 24-1.0343065-turtle.txt
- N-Triples: 24-1.0343065-rdf-ntriples.txt
- Original Record: 24-1.0343065-source.json
- Full Text
- 24-1.0343065-fulltext.txt
- Citation
- 24-1.0343065.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.24.1-0343065/manifest