S I M U L A T I O N O F E L E C T R O M A G N E T I C T R A N S I E N T S IN U N D E R G R O U N D C A B L E S W I T H F R E Q U E N C Y - D E P E N D E N T M O D A L T R A N S F O R M A T I O N M A T R I C E S By LUIS MARTI Ing. Elec. Central University of Venezuela, 1979. M.A.Sc. The University of British Columbia, 1982. A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA November 1986 © Luis Marti, 1986 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia 1956 Main Mall Vancouver, Canada Department V6T 1Y3 Date / j a r r ^ s r /ft? DE-6(3/81) ii A B S T R A C T This thesis presents a new model to simulate the behaviour of underground cable systems under transient conditions. The new cable model belongs to the class of time-domain, frequency-dependent models, and it is directly compatible with the solution algorithm of the EMTP (Electromagnetic Transients Program). The most important feature of the new model is that it takes into account the frequency dependence of the modal transformation matrices and cable parameters, thus overcoming the main limitation of currently-used transmission line and cable models, which assume that the modal transformation matrices are constant Conceptually, the new model is relatively simple. The system parameters which define the behaviour of an underground cable (namely the modal characteristic admittance matrix, the modal propagation matrix, and the modal transformation matrix), are expressed in closed form by approximating them with rational functions in the frequency domain. Therefore, in the time domain, all numerical convolutions can be expressed recursively. The host transients program (to which the model is interfaced) sees the new model as a constant, real admittance matrix, in parallel with a continuously-updated vector current source. The accurate approximation by rational functions of the modal transformation matrix is possible when its elements are continuous and smooth functions of frequency. Standard eigenvalue/eigenvector algorithms are not well suited for this purpose. Therefore, a new procedure to generate eigenvalues and eigenvectors has been developed. This procedure is based on the Jacobi method, and it produces the desired smooth functions of frequency. iii This manuscript presents a number of simulations where the performance of the new cable model is compared with exact analytical solutions. These simulations show an excellent agreement between analytical and numerical answers. The effects of not taking into account the frequency dependence of the modal transformation matrices is illustrated with the simulation of a line-to-ground fault on a three-phase cable. The response of the new cable model is also compared with results measured in a field test The new cable model is numerically stable. Its computational speed is comparable to that of frequency-dependent line models with constant transformation matrices. The new cable model is general. Its extension to the simulation of multiple-circuit overhead transmission lines should also be of considerable practical importance. iv Table of Contents Abstract ii Table of Contents iv List of Tables vi List of Figures vii Dedication xii Acknowledgements xiii Nomenclature 1 INTRODUCTION 2 CHAPTER 1: Mathematical Considerations 7 1.1 Description of the New Cable Model 7 1.2 Synthesis of the Transfer Functions 17 1.3 Synthesis of the Modal Transformation Matrix 21 1.4 Synthesis of the Modal Characteristic Admittance Matrix 39 1.5 Synthesis of the Modal Propagation Matrix 46 1.6 Diagonalization of YZ 57 CHAPTER 2: Numerical Results 60 2.1 Open and Short-Circuit Responses in the Frequency Domain 61 a) Symmetric Excitation 64 b) Zero-Sequence Excitation 74 2.2 Impulse and Square-wave Responses in the Time Domain 83 a) Impulse Response 84 b) Square-Wave Response 89 2.3 Simulation of Single-Phase Line-to-Ground Fault 93 2.4 Comparison with a Field Test 99 2.5 Numerical Stability 104 2.6 Computational Speed 106 CHAPTER 3: Future Research in Cable Modelling and Related Areas 108 3.1 Simulation of Pipe-Type Underground Cables 108 3.2 Extension of the New Cable Model to Overhead Transmission Lines 108 CONCLUSIONS 112 APPENDIX I: Evaluation of Cable Parameters 114 LA Shunt Admittance and Series Impedance Matrices in Loop Quantities 114 LB Transformation Between Loop and Phase Quantities 118 I.C The Generalized Eigenproblem 120 V APPENDIX II: Numerical Recursive Convolutions 123 APPENDIX III: Physical Data of Test Cables 125 References 127 vi List of Tables Table 1: Order of the synthesized elements of Q 30 Table 2: Symmetry within the elements of Q 30 Table 3: Order of the synthesized elements of Y' and A'. Cable length = 10 km c. 39 Table 4: Relative CPU times. EMTP time-step loop 107 vii List of Figures Fig. 1: Multiphase underground cable system 7 Fig. 2: Voltages and currents in an n-conductor transmission system 8 Fig. 3: Single-filament representation of an n- conductor transmission system 9 Fig. 4: Equivalent circuit in the time domain 14 Fig. 5: Magnitude of the elements of column 5 of Q. Evaluated using DCEIGN 23 Fig. 6: Magnitude of the elements of column 6 of Q. Evaluated using DCEIGN 23 Fig. 7: Magnitude of the elements of column 1 of Q. Evaluated using seeded JACOBI 32 Fig. 8: Relative error function: Magnitude, column 1 of Q. Evaluated using seeded JACOBI 32 Fig. 9: Phase angle (in degrees) of the elements of column 1 of Q. Evaluated using seeded JACOBI 33 Fig. 10: Absolute error function: Phase angle (in degrees), column 1 of Q Evaluated using seeded JACOBI 33 Fig. 11: Magnitude of the elements of column 2 of Q. Evaluated using seeded JACOBI 34 Fig. 12: Phase angle (in degrees) of the elements of column 2 of Q. Evaluated using seeded JACOBI 34 Fig. 13: Magnitude of the elements of column 3 of Q. Evaluated using seeded JACOBI 35 Fig. 14: Phase angle (in degrees) of the elements of column 3 of Q. Evaluated using seeded JACOBI 35 Fig. 15: Magnitude of the elements of column 4 of Q. Evaluated using seeded JACOBI 36 Fig. 16: Phase angle (in degrees) of the elements of column 4 of Q. Evaluated using seeded JACOBI 36 Fig. 17: Magnitude of the elements of column 5 of Q. Evaluated using seeded JACOBI 37 Fig. 18: Phase angle (in degrees) of the elements of column 5 of Q. Evaluated using seeded JACOBI 37 viii Fig. 19: Magnitude of the elements of column 6 of Q. Evaluated using seeded JACOBI 38 Fig. 20: Phase angle (in degrees) of the elements of column 6 of Q. Evaluated using seeded JACOBI 38 Fig. 21: Magnitude (in S) of element 2 of Y ' £ 41 Fig. 22: Relative error function: Magnitude of element 2 of Y' 41 Fig. 23: Phase angle of (in degrees) of element 2 of Y'c 42 Fig. 24: Absolute error function: Phase angle (in degrees) of element 2 of Y ' c . . . . 42 Fig. 25: Magnitude (in S) of elements 1 and 2 of Y' 43 Fig. 26: Phase angle (in degrees) of elements 1 and 2 of Y ' c 43 Fig. 27: Magnitude (in S) of elements 3 and 4 of 44 Fig. 28: Phase angle (in degrees) of elements 3 and 4 of Y ' c 44 Fig. 29: Magnitude (in S) of elements 5 and 6 of Y ' £ 45 Fig. 30: Phase angle (in degrees) of elements 5 and 6 of Y ' c 45 Fig. 31: Phase angle (in degrees) of element 6 of A' 47 Fig. 32: Magnitude of element 6 of A' 47 Fig. 33: Relative error function: Magnitude of element 6 of A' 48 Fig. 34: Phase angle (in degrees) of element 6 of A' after extracting COT 48 Fig. 35: Absolute error function: Phase angle (in degrees) of element 6 of A' 49 Fig. 36: Magnitude of elements 1 and 2 of A' 53 Fig. 37: Phase angle (in degrees) of elements 1 and 2 of A' after extracting COT 53 Fig. 38: Magnitude of elements 3 and 4 of A' 54 Fig. 39: Phase angle (in degrees) of elements 3 and 4 of A' after extracting COT 54 Fig. 40: Magnitude of elements 5 and 6 of A' 55 Fig. 41: Phase angle (in degrees) of elements 5 and 6 of A' after extracting COT 55 Fig. 42: Magnitude of element (2,1) of A (phase) 56 ix Fig. 43: Magnitude of element (2,2) of A (phase) 56 Fig. 44: Single-phase transmission system 61 Fig. 45: Three-phase underground cable. Symmetric excitation 64 Fig. 46: Open-circuit response (magnitude). Symmetric excitation: Core 1 (new model) 66 Fig. 47: Open-circuit response (magnitude). Symmetric excitation: Core 1 (constant Q, evaluated at 0.1 Hz) 67 Fig. 48: Open-circuit response (magnitude). Symmetric excitation: Core 1 (constant Q, evaluated at 60 Hz) 67 Fig. 49: Open-circuit response (magnitude). Symmetric excitation: Core 1 (constant Q, evaluated at 5 kHz) 68 Fig. 50: Short-circuit response (magnitude). Symmetric excitation: Core 1 (new model) 69 Fig. 51: Short-circuit response (magnitude). Symmetric excitation: Core 1 (constant Q, evaluated at 0.1 Hz) 70 Fig. 52: Short-circuit response (magnitude). Symmetric excitation: Core 1 (constant Q, evaluated at 60 Hz) 71 Fig. 53: Short-circuit response (magnitude). Symmetric excitation: Core 1 (constant Q, evaluated at 5 kHz) 72 Fig. 54: Short-circuit response (phase angle). Symmetric excitation: Core 1 (new model) 73 Fig. 55: Short-circuit response (phase angle). Symmetric excitation: Core 1 (constant Q, evaluated at 5 kHz) 73 Fig. 56: Three-phase underground cable. Zero-sequence excitation 74 Fig. 57: Open-circuit response (magnitude). Zero-sequence excitation: Core 1 (new model) 76 Fig. 58: Open-circuit response (magnitude). Zero-sequence excitation: Core 1 (constant Q, evaluated at 0.1 Hz) 77 Fig. 59: Open-circuit response (magnitude). Zero-sequence excitation: Core 1 (constant Q, evaluated at 60 Hz) 77 Fig. 60: Open-circuit response (magnitude). Zero-sequence excitation: Core 1 (constant Q, evaluated at 5 kHz) 78 Fig. 61: Short-circuit response (magnitude). Zero-sequence excitation: Core 1 (new model) 79 X Fig. 62: Short-circuit response (magnitude). Zero-sequence excitation: Core 1 (constant Q, evaluated at 0.1 Hz) 80 Fig. 63: Short-circuit response (magnitude). Zero-sequence excitation: Core 1 (constant Q, evaluated at 60 Hz) 81 Fig. 64: Short-circuit response (magnitude). Zero-sequence excitation: Core 1 (constant Q, evaluated at 5 kHz) 82 Fig. 65: Three-phase, 230 kV underground cable 83 Fig. 66: Fourier series approximation of a voltage impulse. N = 500, T=lms 86 Fig. 67: Impulse response, receiving-end voltage. Core 1 87 Fig. 68: Impulse response, receiving-end voltage (third cycle). Core 1 88 Fig. 69: Impulse response, receiving-end voltage (third cycle). Sheath 1 88 Fig. 70: Fourier series approximation of a square wave. N = 1000, T=40ms 90 Fig. 71: Square-wave response, receiving-end voltage. Core 1 91 Fig. 72: Square-wave response, receiving-end voltage (second cycle). Core 1 92 Fig. 73: Square-wave response, receiving-end voltage (second cycle). Sheath 1 92 Fig. 74: Simulation of a single-phase line-to-ground fault 93 Fig. 75: Receiving-end voltage. Core 3 (sheaths are grounded) 95 Fig. 76: Receiving-end voltage. Core 2 (sheaths are grounded) 95 Fig. 77: Receiving-end voltage. Core 1 (sheaths are grounded) 96 Fig. 78: Receiving-end current Core 3 (sheaths are grounded) 96 Fig. 79: Receiving-end voltage. Core 3 (sheaths are not grounded) 97 Fig. 80: Receiving-end voltage. Core 2 (sheaths are not grounded) 97 Fig. 81: Receiving-end voltage. Core 1 (sheaths are not grounded) 98 Fig. 82: Receiving-end current Core 3 (sheaths are not grounded) 98 Fig. 83: Field test connection diagram 99 Fig. 84: Sending-end voltage. Core 2 (input waveshape 0 X40MS) 101 Fig. 85: Second crossbonding- point voltage. Sheath 1 (input waveshape 0 X40MS) . . . . 101 xi Fig. 86: Second crossbonding-point voltage. Sheath 2 (input waveshape OxAOus) . . . . 102 Fig. 87: Second crossbonding-point voltage. Sheath 3 (input waveshape 0 X40MS) . . . . 102 Fig. 88: Sending-end voltage. Core 2 (input waveshape 0X45MS) 103 Fig. 89: Magnitude of column 1 of Q of a double-circuit transmission line. Two 500 kV circuits sharing the same tower Ill Fig. 90: Magnitude of column 1 of Q of a double-circuit transmission line. 115 kV and 230 kV circuits sharing the same tower Ill Fig. 1.1: Cross-section of an underground cable ' 115 Fig. 1.2: Single-phase underground cable 116 Fig. III.l: Cable 1 configuration 126 Fig. III.2: Cable 2 configuration 126 xii To Chris xiii Acknowledgements I would like to express my appreciation to all the people who contributed, directly or indirectly, to the realization of this thesis project. In particular I would like to thank Professor H. W. Dommel, for his encouragement, advice, and patience; Professor J. R. Marti, who provided and modified FI.T-S to suit the purposes of this thesis, for his encouragement, advice, and many stimulating discussions; Professor A. Ametani, for providing the cable data used in the comparison with the field test; Chris, Jose, and Mae, for revising this manuscript. Also, the financial support of British Columbia Hydro and Power Authority, the University of British Columbia, and Bonneville Power Administration, is gratefully acknowledged. 1 Nomenclature The following conventions will be followed in this thesis: - Upper case letters indicate "frequency-domain" quantities, whereas lower case letters are used for their "time-domain" counterparts (e.g., V in the frequency domain and v(t) in the time domain). Unless otherwise indicated, the Fourier Transform is assumed for the transformation between time and frequency domains. - Primed symbols are used to refer to "modal" quantities; otherwise, "phase" components should be assumed (e.g., v'(t) for modal voltage, and v(t) for phase voltage). - All equations are assumed to be matrix equations. Equations describing frequency-domain relationships contain square, complex matrices dimensioned nxn, and complex vectors of length n. Equations describing time-domain relationships only contain square, real matrices dimensioned nxn, with real vectors of length n. - Subscript "k" is used to indicate quantities measured at the sending end, while subscript "m" is used for receiving-end quantities. - The term "transmission line" will be used to refer to high-voltage overhead transmission lines. The term "cable" will be used to refer to high-voltage buried conductors consisting of a core and one or more concentric sheaths. Whenever necessary, deviations from these conventions will be indicated explicitly. 2 INTRODUCTION The accurate modelling of underground high-voltage cables becomes increasingly important as the load concentration in populated areas increases, and buried cables often become the only viable form of power delivery for aesthetic, safety, and economic reasons. Insulation failures caused by transient overvoltages on transmission lines are usually self-restoring. After a flashover across an insulator, the arc will be extinguished when the line is temporarily switched off, and normal operation can be resumed after re-closing. In contrast, a fault in an underground cable can result in irreversible damage to its insulation, which is not self-restoring. Therefore, the affected section must be located and replaced before the cable can be put back into service. To reach a suitable compromise between high insulation levels at high costs and lower insulation levels at lower costs (but potentially more expensive if failures occur), the behaviour of a cable system under abnormal operating conditions must be predicted. Many of the analytical tools presently used in the simulation of transient phenomena in underground cables are based on models originally developed for transmission lines. However, many of these models are not accurate enough for underground cables, because some of the simplifying assumptions in the line models are not valid for cables. The existing concern in the accurate modelling of underground cables is reflected by the numerous models proposed in the technical literature (e.g., [1] to [6], to mention a few). The following is a brief overview of models which, to the author's knowledge, are the most widely used in the simulation of transient phenomena in underground cable systems. 3 Overview Before the widespread use of digital computers, specialized analog computers in the form of "Transient Network Analyzers" (TNA) were used for the simulation of transient phenomena in power systems. In a TNA study, a cable is represented with a number of cascaded it-sections whose parameters (R and L) are evaluated at a single frequency. A more sophisticated form of this model represents the ground return impedance with a suitable combination of two or more parallel R-L branches [8]. The validity of this type of lumped-parameter representation is restricted to relatively short cables, and, in general, its frequency response is only good in the vicinity of the frequency at which the parameters are evaluated. The increasing availability of digital computers has resulted (over the past decade or so) in a shift from analog to digital simulations. The accuracy and flexibility of digital programs have made them the dominant tool in modern transient simulations. Two basic types of (digital) cable models can be distinguished, according to the solution techniques used by their host programs: 1) Time-domain models 2) Frequency-domain models. 1) Time-domain models In this class of models the solution of the differential equations which describe the transient behaviour of a transmission system is carried out in the time domain without the explicit use of inverse Fourier or Laplace transforms. Over the past two decades, a relatively large number of time-domain transmission line models have been proposed. Among those which have been used in the simulation 4 of underground cables, two models deserve special mention: 1) Lumped-parameter models. Conceptually, this type of model is very similar to the TNA representation; that is, the cable is represented by cascaded ir-sections evaluated at a given frequency [9]. Although lumped-parameter models are computationally fast, they do not take into account the distributed nature of the parameters, and their frequency response is generally restricted to the neighbourhood of the frequency at which R and L are evaluated. ii) Distributed-parameter, frequency-dependent models. In this type of model, the differential equations are decoupled using modal transformation matrices. The frequency dependence and the distributed nature of the cable parameters are taken into account [5]. However, the validity of this model is restricted because the modal transformation matrices are assumed to be real and constant This assumption, although valid in the case of most overhead transmission lines, leads to poor results in the case of underground cables because their modal transformation matrices depend very strongly on frequency [30]. 2) Frequency-domain models. In this class of models, the response of the transmission system is evaluated direcdy in the frequency domain. The time-domain solution is found using inverse transformation algorithms such as the EFT (Fast Fourier Transform) [6]. The frequency dependence of the cable parameters and of the modal transformation matrices is taken into account Even though inherent numerical problems such as aliasing and Gibbs' oscillations have been alleviated for most practical purposes (using windows and other 5 oscillation-suppressing techniques), the applicability of these models tends to be restricted by the limitations of their host programs. To date, there are no general-purpose transient analysis programs in the frequency domain with the overall simulation capabilities of time-domain programs such as the EMTP [10]. Sudden changes in the network configuration (such as faults, opening and closing of circuit breakers, etc.) and the modelling of non-linear elements cannot be handled easily with frequency-domain solution methods. It is then the limitation of the host program, rather than the model's deficiency, what has prevented the widespread use of frequency-domain cable models; nevertheless, these models are often used as benchmark solutions to verify the validity of time-domain cable models, and in the simulation of relatively simple network configurations. Also, frequency-domain models are very efficient in the simulation of large numbers of crossbonded cable sections, as long as no additional network components are connected at the end points of each section. The new cable model The cable model presented in this thesis belongs to the class of time-domain, frequency-dependent models. It overcomes the main limitation of existing time-domain line models; that is, it takes into account the frequency dependence of the modal transformation matrices.1 Also, by being compatible with the solution algorithm of the EMTP, it enjoys all the implicit advantages of a versatile host program. Therefore, it also overcomes the main limitation of frequency-domain cable models. 1 A time-domain cable model which takes into account the frequency dependence of the transformation matrices has also been developed at Doshisha University in Japan by H. Nakanishi and A. Ametani, in parallel with the work reported in this thesis. Their recent publication [7] indicates that the two approaches are different in many aspects. A careful investigation of these two models would be a worthwhile project for future work on cable modelling. 6 The new cable model is accurate, as demonstrated with comparisons with analytical and field test results. From a computational point of view, the speed of this model is comparable to the speed of frequency-dependent models with constant transformation matrices. Theoretically, the new cable model is applicable to any transmission system. However, the research effort of this thesis project has only been focused on the simulation of underground cables. The future extension to the simulation of multiple-circuit overhead lines should also be of considerable practical importance. 7 CHAPTER 1 Mathematical Considerations 1.1 Description of the New Cable Model Consider a multiphase, underground cable system of length / such as the one illustrated in Figure 1. Node k Node m ^1 ] ^ able #1 y}j J ) Cable #2 ) Cable #j 1 — c > 0 I x Fig. 1: Multiphase underground cable system. This transmission system contains a total of "n" conductors, including sheaths and armours, if present Voltages are measured with respect to ground; sending-end currents flow into the cable, while receiving-end currents out of the cable, as shown in Figure 2. A more compact representation is shown in Figure 3, where V and I represent voltage and current vectors of dimension n. 8 In the frequency domain, the behaviour of this n-conductor system can be described with the well-known differential equations dV x dx dl X dx = ZI = YV (1) (2) where V and I are the voltages and currents measured at a distance x from the sending end; Z and Y are the series impedance and shunt admittance matrices, respectively. Note that Z and Y are nxn, complex, symmetric matrices, while and I are complex vectors of dimension n. Node k Node m //////////////////////////////////////////////////////////////////// Fig. 2: Voltages and currents in an n-conductor transmission system. 9 * o-A -o — O m ///////////////////////////////////////////////////// Fig. 3: Single-filament representation of an n-conductor transmission system. The solution of equations (1) and (2), in terms of sending and receiving-end quantities, can be expressed as Y„ V +1 = F = A F, c m m m k (3) Y c V * " h = B * = A Bm- (4) where Y c = »^YZT-1Y (5) A = exp{-/YZ"/l. (6) Y £ is the characteristic admittance matrix, and A is defined as the propagation matrix; F and B are intermediate vector functions. In the time domain, F and B can be visualized as waves travelling in forward and backward directions, respectively. Let us assume that YZ is diagonalizable throughout the entire frequency range of interest to transient calculations (typically from dc to 1 MHz). Then, there exists a non-singular matrix Q, such that 10 Q - ' Y Z Q = D^, where is a diagonal matrix whose non-zero entries are the eigenvalues of YZ, and the columns of Q are its associated eigenvectors (see Appendix I.C). From equation (6), it is then clear that the propagation matrix A can be diagonalized or transformed into the modal domain; that is, A' = Q" 1 A Q, (7) where A' (the modal propagation matrix) is a diagonal matrix given by A'= exp{-/JT^ /}. yz Note that the diagonal entries of A' are the eigenvalues of A. The characteristic admittance matrix can also be transformed into the modal domain (i.e., expressed in diagonal form) with the relationship Y' c = Q " 1 Y c Q - T , (8) where Q" T denotes the transpose of Q" 1 . Note that, in general, the diagonal entries of Y'c are not the eigenvalues of Y . Introducing equations (7) and (8) into (3) and (4) we obtain, after some algebraic manipulation, Y' V' + I' = F ' = A'F'. (9) c m m m k where 11 V = Q T V r = Q" 1 i ( i i ) F' = Q- 1 F B' = Q" 1 B. Equations (9) and (10) represent a system of n decoupled linear equations, where Y' and A' are complex-valued functions of frequency. Transformation of equations (9) and (10) into the time domain gives y'c(t) • v^ (t) - i^ (t) = b'^ (t) = a'(t) • b'Jt). (13) where the relationship between modal and phase voltages is now given by v'(t) = qT(t) • v(t) (14) i'(t) = q- 1(0 * i(t) (15) and the symbol "*" is used to indicate matrix-vector convolutions (a matrix-vector convolution can be visualized as a standard matrix-vector product where all multiplications are replaced by scalar convolutions). Consider now equation (12). Let us assume that the elements of y^ (t) can be expressed as a finite sum of exponentials of the form 12 m p(t) = k06(t) + I k exp{-pt} u(t), 1=7 (16) where 5(t) is the Dirac impulse function; u(t) is the unit step function; k0, k^ . and p^ . are real constants. The numerical convolution of p(t) with a given function v(t) can then be evaluated using well-known recursive techniques [1]. Therefore, as indicated in Appendix II, y'(t) * v' (t) in equation (12) can be expressed as where y' is a real, diagonal matrix, and h' (t) is a function of past history values of CO 171 y g^ (t) and v (^t) (that is, values already known prior to time t). Assuming that the elements of the modal transformation matrix (in the time domain) can also be expressed in the exponential form shown in equation (16), we can evaluate the convolution indicated in equation (14) for the receiving-end voltages; that is, y'(t) * v'(t) = g'(t) = y' v'(t) + h' (t), (17) v'(t) = q 0 T v (t) + h' (t). (18) Similarly, from equation (15) we can evaluate the receiving-end currents, q 0i'(t) + h' (t) M U m v ' m 3 v ' - *>-1h;3(t), (19) 13 where q 0 is a real, nxn matrix; h' (t) is a function of past history values of v (t) " I 2 rn and v'(t); and h' (t) is a function of past history values of i (t) and i'(t). rrv' rrv' rrv' The only convolution left to evaluate in equation (12) involves the modal propagation function a'(t). Let us assume that a'(t) can be expressed as a Finite sum of exponentials of the form where k;, p;. and r are real constants. Note that the time-displacement constant r is numerically close to the travel time of the fastest frequency component of a wave propagating on a given mode. It is shown in Appendix II that the numerical convolution of p(t) with a given time function f(f) depends exclusively on past history. Therefore, f^ (t) = a'(t) * f^ (t) in equation (12) is completely defined, at time t, in terms of pre-recorded past history values of f^ (t) and f^ (t). Introducing equations (17) to (19) into (12), we then obtain m p(t) = I k. exp{-p/t-T)l u(t-T), (20) (21) where = 4° y' C 0 T h (t) = q„[f'(t) - y' h' (t) - h' (t)] + h' (t). (22) Following an analogous procedure for equation (13), we obtain (23) 14 where tyt) = Qo [b^ (t) - y' c o h'k (t) - h'k (t)] - h'^ (t). (24) Note that yg^ is a real, symmetric (nxn) matrix; and that h (^t) and h t^) are defined as equivalent history current sources. At time t, these vectors are completely defined in terms of variables known from previous time steps. Equations (21) and (23) can be represented by the equivalent circuit shown in Figure 4. k O 0 vfc(t) eq k i (0 m h m ( t ) % -O m v (t) m Fig. 4: Equivalent circuit in the time domain. This equivalent circuit is compatible with the solution algorithm of the EMTP. In fact, the EMTP representation of lossless overhead lines [10] and the representation of frequency-dependent overhead lines with constant transformation matrices [5] also share the same form; that is, a real, constant admittance matrix in parallel with an equivalent vector current source. From a computational point of view, only the updating of these equivalent sources is different; this makes the implementation of the interface between the EMTP and the new cable model relatively straightforward and completely transparent as far as the main core of the host program is concerned. 15 The equivalent circuit for frequency-dependent overhead line models with constant transformation matrices can be derived as a particular case of the general underground cable model presented here. If q(t) is assumed to be constant (e.g., q(t) = q0), then h' (t), h', (t), h' (t) and h', (t) will become zero (see equations (18) and (19)), and hw(t) = Qo[r(t) - h; i (t)] (25) "eq Vk® ~ W = h*« h^O = q0[b'^ (t) - h^(t)]. (26) This result is important from a computational point of view, because it permits the implementation of only one set of subroutines to model overhead transmission lines (if the use of a constant Q is desirable) as well as underground cables. In summary, the derivation of the new cable model is conceptually simple: - The solution of the differential equations which describe an n-conductor transmission system in the frequency domain (equations (3) and (4)) is transformed into the modal domain (equations (9) to (11)) in order to obtain a decoupled system of linear equations. - The time-domain solution is then formulated through direct use of the inverse Fourier Transform. What used to be products in the frequency domain will now become convolutions in the time domain (equations (12) to (15)). - All numerical convolutions between system variables (e.g., v(t), i(t), etc.) and transfer functions (a'(t), y^ ,(t) and q(t)) are expressed recursively, given that all transfer functions, including the modal transformation matrix q(t), can be expressed as finite 16 sums of exponentials (equations (16) and (20)). - Algebraic manipulation eventually reduces the equations to a form which is directly compatible with the solution algorithm of the EMTP (equations (21) and (23)). Clearly, the most important assumption in the derivation of the new model is that y^ .(t), a'(t) and q(t) can be expressed in closed form. Therefore, the justification for this assumption will be the main concern of the remaining sections of this chapter. 17 1.2 Synthesis of the Transfer Functions In the derivation of the equivalent circuit of the new cable model, it was assumed that the elements of matrices a'(t), q(t), and y'c(i) can be expressed as finite sums of exponentials. In the case of q(t) and y^ ,(0, this exponential representation has the form m p(t) = k06(t) + Z k exp{ -p t} u(t), i=l while in the case of a'(t), its elements are expressed as m p(t) = Z k exp{-p/t-r)l u(t-r). 1=1 In the frequency domain, the series representation of the elements of y^ (t) and q(t) has the form m k. P(CJ) = k 0+ Z -i=l jcj + p;. while for the elements of a'(t) R(CJ) = P(co) exp{-jo)T} m k. P(w) = Z l -If the physical characteristics of an underground cable are available, there are well-known formulae which permit the calculation of the Z and Y matrices at a given frequency (see Appendix LA). 18 Once Z and Y are known, the matrix product YZ. must be diagonalized in order to find Q, Y' c and A'; that is, the following eigenvalue/eigenvector problem must be solved: YZ Q = Q from which A' = exp{-/TJT /} y2. Y' = / ( D l " 1 D : c v yr y where D = Q - ' Y Q ' T . y D is a diagonal matrix containing the eigenvalues of YZ; the columns of Q are the corresponding eigenvectors, and matrices Y'c, A' and are diagonal. When this process is carried out at a sufficiently large number of frequencies, the elements of Y', A' and Q can be obtained as smooth functions of frequency evaluated at discrete points. Typically, ten points per decade, starting at a frequency of 0.001 Hz, is adequate. The transition between the discretized functions of frequency and the representation by rational functions required by the new cable model is made with a curve-fitting algorithm which will be referred to as FIT-S. This curve-fitting algorithm was developed by J. R. Marti in 1981 to produce approximations by rational functions of the transfer functions of overhead transmission lines with constant transformation matrices. With some (albeit minor) contributions on the part of the author, the algorithm has been refined and updated by J. R. Marti in order to satisfy some of the particular needs of the 19 functions encountered in this thesis project Given a complex-valued function of frequency, ? e x a c t , FIT-S traces its magnitude function with the magnitude function of P f l^- 0 J C. where m k. p » p = k 0+ I (27) exact approx " . . . v ' i=l }co + p / FIT-S forces the poles and zeros of the approximating function to be real and negative; therefore, P belongs to the class of minimum-phase-shift functions. This approx r implies that if the magnitude function |P f l is known, then its phase function arglP^^^Hs uniquely defined [11]. Therefore, if ^ e x a c t is a minimum-phase-shift function, the condition 1 exact 1 approx is sufficient to guarantee that the phase functions will also match; that is, arg{P J arg{P }. ° exact e approx The accuracy of FIT-S depends on the number of terms used in P . and r appro* on the shape of the magnitude of the function to be approximated. For example, to obtain errors under 4% in the approximation of the significant magnitudes of the elements of Y', A' and Q for the three-phase, 230 kV underground cable described in Appendix III, the average number of terms required was 18, 17 and 7, respectively1 Note that the order of the approximating function can be controlled to a certain extent FIT-S gives the best approximation within a specified number of terms. For the 'This underground cable will be used in all the numerical examples shown in this thesis, unless explicitly indicated otherwise. 20 example mentioned above, the maximum number of poles permitted was 25. Also note that the errors in arg { P } are of the same order of magnitude as the errors in approx ° IP J . as long as P . is a minimum-phase-shift function. 1 approx exact From the above discussion, it is clear that to obtain a good match in both magnitude and phase angle, the functions to be approximated must be minimum-phase-shift In the case of the underground cables studied in this thesis project, the elements of Y', A' and Q can be manipulated into minimum-phase-shift form. A detailed discussion of the procedures and difficulties involved is given in the next sections. 21 1.3 Synthesis of the Modal Transformation Matrix In order to obtain the eigenvalues of YZ, it is standard practice to solve the eigenproblem Note that although Z and Y are symmetric matrices, their product is not symmetric. There are several numerical algorithms which can evaluate the eigenvalues and eigenvectors of non-symmetric, complex matrices. One of the most widely used algorithms is available in subroutine form under the name of DCEIGN.2 In the past, DCEIGN has been used in EMTP support routines such as LINE CONSTANTS [13]. However, DCEIGN is not well suited for the evaluation of the eigenvectors of YZ as smooth functions of frequency. Its main drawbacks are: a) When the eigenproblem is solved at a given frequency, DCEIGN cannot "remember" the relative order of the eigenvectors obtained in the previous frequency step. For example, suppose that the transformation matrices obtained at two consecutive frequency YZ Q = Q D steps are 1.00 0.00 Q("o) 0.00 1.00 1.20 0.10 Q("i) 0.10 0.98 DCEIGN, however, could have produced the following matrix at u> = CJ , 2The solution algorithm used in DCEIGN is based on the Householder transformation and subsequent QR decomposition. This subroutine is part of the mathematical library package EISPACK [12]. 22 -0.10 1.20 Q > , ) = 0.98 0.10 where the relative positions of the columns of Q'(o),) have been switched. It could be argued that a sorting algorithm could detect the switchover; however, it is not possible to obtain consistent scaling factors for the columns of Q, unless their relative positions are known. For example, DCEIGN could also have produced the following matrix A tracking algorithm could decide that an eigenvector switchover has occurred, when in reality, only the scaling factor of column 1 of Q"(o>1) has changed. In a more realistic case, the problem is considerably more complicated, because there seldom are only two columns to choose from, and because the elements of Q are complex numbers. b) When the eigenvalues of YZ are very closely packed, DCEIGN shows numerical stability problems. If, for example, Y and Z are balanced,3 then YZ is also balanced and DCEIGN collapses [13]. Figures 5 and 6 show the magnitude of the elements of columns 5 and 6 of Q evaluated with DCEIGN. Note the numerous eigenvector crossovers up to 100 kHz. In the 100 kHz to 10 MHz region, crossovers are combined with numerical instability due to closely packed eigenvalues. -1.20 -0.10 1 Q"(c_.,) = -0.10 0.98 3 A matrix is said to be balanced when all its diagonal elements are equal among themselves, and all its off-diagonal elements are also equal among themselves. Such a matrix would only have two distinct elements. 1 . 2 - , F r e q u e n c y (Hz) Fig. 5: Magnitude of the elements of column 5 of Q. Evaluated using DCEIGN. F r e q u e n c y (Hz) Fig. 6: Magnitude of the elements of column 6 of Q. Evaluated using DCEIGN. 24 Clearly, the elements of Q cannot be approximated as shown. The possibility of a tracking algorithm was explored, and to a large extent it was possible to unscramble the eigenvectors. However, the reliability was precarious, at best, and the instability problem at high frequencies remained unresolved. Fortunately, the difficulties outlined above have been solved by reformulating the eigenproblem YZ Q = Q in the form where H R = R D^; (28) R = ( / ( D ^ T 1 S"1 Q (30) D . = S" 1 Y S. (31) In other words, instead of finding the eigenvectors of YZ, we find the eigenvectors of a symmetric matrix H which has the same eigenvalues as YZ [14]. The eigenvector matrix R (associated with the new matrix H), on the other hand, is related to the eigenvector matrix of the original problem by equation (30). Note that D is a yeig diagonal matrix containing the eigenvalues of Y, and S is its corresponding eigenvector matrix (see Appendix I.C). 25 The solution of the modified eigenproblem has the following advantages: a) Since matrix H is symmetric, its eigenvalues and eigenvectors can be evaluated using an extension of the Jacobi algorithm for complex symmetric matrices. The Jacobi algorithm is a completely stable method for finding the eigenvalues and eigenvectors of a real symmetric matrix [15]. In this thesis project, this algorithm was modified to solve the eigenproblem for the complex symmetric case, and implemented in subroutine form under the name JACOBI. The use of the Jacobi method for the complex symmetric case does not appear to be considered in the literature; therefore, a formal proof of its stability in the complex case could not be found. However, in the cases studied in this thesis, JACOBI always converged to the correct solution even in the cases where DCEIGN failed. The modified Jacobi algorithm finds the eigenvalues and eigenvectors of a symmetric matrix by successive rotations of its would-be eigenvector matrix. If a good initial guess of the eigenvector matrix R is known, convergence to the true eigenvector matrix is very fast because the extent of the rotations is small. Another way to visualize this, is that if the matrix to be diagonalized is already quasi-diagonal (because it has already been pre and post-multiplied by a good guess of its eigenvector matrix), disturbances or rotations acting on the initial guess are minimal. This property is ideally suited to our particular needs because a good estimate of the eigenvector matrix is already known: the eigenvector matrix obtained in the previous frequency step. For example, assume that R(CJ 0) is known; that is, R - 1 (CJ 0 ) H(o)0) R(w0) = D > 0 ) ; (32) 26 and let us also assume that we wish to diagonalize H(CJ ,), where OJ , is sufficiently close to w 0. 4 Since the variation of the elements of H between a> = a>0 and o) = a>, is smooth, we have that R" 1 (o>o) H(_>,) R(CJ 0) = H , , (33) where H, is practically diagonal. If we use JACOBI to find the eigenvalues and eigenvectors of H, R i 1 H, R, = D^Cw,). we will find that R, is very close to the identity matrix and that the convergence to the true solution is very fast Once the eigenvectors of H, are found, the eigenvectors of H(a>, ) can be obtained using the relationship R(_),) = R, R(_>0). This "seeding" process is computationally fast, despite the seemingly large number of matrix products involved. Furthermore, the eigenvector switchover problem is altogether solved without the need for a sorting algorithm. Note, however, that seeding alone is not sufficient to prevent switchovers: if DCEIGN is used with seeding, or if JACOBI is used without seeding, eigenvector switchover occurs. b) The eigenvector matrix R satisfies ST S = D [17]. In particular, R can be scaled so that D = U, where U is the identity matrix. This results into two side benefits: i) When H{u^) is seeded to find H,, R(a>0) does not have to be inverted because R (CJ 0 ) - 1 =R (OJ 0) T. 4When R(u>0) is not known because <_>, is the first frequency at which R(CJ) is calculated, then R(CJ 0) is assumed to be the unit matrix. 27 ii) If the variation of the elements of R from a>0 to u, is smooth, then R(w0)R( u ( ) T =U", where U'is very close to the identity matrix as long as no eigenvector switchover occurs between C J 0 and CJ^. In the case of the underground cables studied, U-U'is practically the zero matrix; that is, the magnitude of the elements of U - U' is generally 0.001 or less, depending on the rate of change of the elements of R in the frequency range in consideration. If, for example, there were a switchover between columns 1 and 3 of R, then the magnitude of elements (1,1), (3,3), (1,3), and (3,1) of U - U' would increase by several orders of magnitude. Therefore, the "complex orthogonality" of R [17] provides the means for a very simple form of switchover checking, given that U-U'will remain small as long as no eigenvector switching occurs. It should be noted that in all cases studied in this thesis, not a single switchover was detected when the modified eigenproblem was solved with JACOBI and seeding. Even when only two points per decade were used to generate Q( CJ) the solution algorithm still produced correct answers. The error-checking provision mentioned above will probably be triggered only in the case of very poorly-behaved eigenvectors. The reformulation of the eigenproblem presented here is general. However, when applied to underground cables of the type shown in Appendix III, some simplifications are possible. If matrices Z and Y are formulated in loop-equation form (see Appendix I.B), is a diagonal matrix; therefore, equations (28) to (31) can be rewritten as follows H, R, = R, D loop loop loop yz 28 RicoP = < ^ y1 Qioop <34> ^loop loop ^ loop loop Note that Y ^ ^ is proportional to CJ, i.e., Y^00^= w(tan8C o^o^ +jC o^o^ ); therefore, it has to be evaluated only once. The transition between loop and node equations is made through a frequency-independent similarity transformation (see Appendix LB); consequently, the eigenvalues of H are the same in both loop and node quantities, while R^^ is related to R through a simple matrix product In the cable constants support routines written as part of this thesis project, all the calculations are made using loop equations. Since A' and Y'c are the same in both loop and node quantities, only has to be transformed into node quantities. In the case of underground cables, no penalty is paid by solving the modified eigenproblem, despite the additional operations needed to seed H. However, computational speed is not a major problem as far as the calculation of eigenvalues and eigenvectors is concerned because it involves less than 5% of the computational effort needed to evaluate cable parameters, which, in turn, represents less than 5% of the total computational effort when the synthesis of rational functions is taken into consideration. Even though measurements of CPU time show that seeded Jacobi is 10% faster than DCEIGN, the overall effect in the evaluation of the cable parameters is practically negligible. Once the matrix function Q(w) has been generated for the entire frequency range of interest, an additional amount of manipulation is necessary before its elements are approximated using FIT-S. The columns of Q are the eigenvectors of YZ; therefore, they can be multiplied by any constant and still be legitimate eigenvectors. If FIT-S is to produce accurate approximations of both the magnitude and the phase angles of the 29 elements of Q, they must be minimum-phase-shift functions. The extra degree of freedom in the selection of the eigenvectors of YZ must then be used to insure that the elements of Q are indeed minimum-phase-shift. In the case of the underground cables studied, it is sufficient to multiply each eigenvector by a factor such that one of its elements becomes constant and real. This scaling process automatically forces the rest of the elements of Q to be minimum-phase-shift functions, and it has the additional advantage that at least n of the elements of Q need not be synthesized with rational functions. In the case of the underground cables studied in this thesis project, it has also been found that there is a considerable amount of symmetry within the elements of Q; that is, many of its elements are identical through the entire frequency range. This reduces the number of approximations to be generated. In the case of our reference cable, only 14 elements had to be synthesized (out of a total of 36, if we consider the entire 6x6 matrix, or out of 30, if we exclude the anchoring elements which were forced to be constant). With the elements of Q properly scaled, we can now proceed to obtain their approximations by rational functions using FIT-S. The results for the above mentioned cable are shown in Figures 7 to 20. Figure 7 shows the magnitude function of the elements of column 1 of Q. Note that this graph contains 12 superimposed curves: the exact values (shown in dashed trace) and their corresponding approximations (shown in solid trace). Also, note the degree of symmetry, as well as the accuracy of the synthesized functions. Figure 9 shows the phase angles for column 1 of Q, and Figures 8 and 10 show the errors in magnitude and phase, respectively. The relative magnitude errors have been calculated using the formula error{%) = 100 [If J - If J]/|f , J v ' 1 exact 1 approxJ 1 exact,max 30 The absolute phase angle errors are calculated as the difference (in degrees) between the phase angles of f „ and f r ° exact approx error(deg.) = arg{ f ,} - arg{ f }. v 6 ' & exact 6 approx Note that the differences between the exact and the synthesized forms of column 1 of Q are larger than the differences observed among the other columns of Q (for this particular cable). col. 1 col. 2 col. 3 col. 4 col. 5 col. 6 row 1 6 9 6 14 0 6 row 2 2 8 11 8 4 0 row 3 0 6 0 10 0 0 row 4 6 0 10 0 0 0 row 5 6 9 6 14 0 6 row 6 2 8 11 8 4 0 1: Order of the synthesized elements of Q. col. 1 col. 2 col. 3 col. 4 col. 5 col. 6 row 1 (1.1) (1.2) (1.3) (1.4) (anchor) (1,6) row 2 (2,1) (2,2) (2,3) (2,4) (2,5) (anchor) row 3 (anchor) (3.2) (anchor) (3,4) (3,5) (3,5) row 4 (4,1) (anchor) (4,3) (anchor) (3,5) (3,5) row 5 (1.1) (1.2) (1.3) (1.4) (anchor) (1.6) row 6 (2,1) (2.2) (2,3) (2,4) (2.5) (anchor) Table 2: Symmetry within the elements of Q. 31 Figures 11 to 20 show the magnitudes and phase angles for the remaining elements of Q, and Table 1 shows the number of poles required in the approximation of each element of Q. Table 2 illustrates the degree of symmetry encountered in this particular case (e.g., elements (5,1) and (1,1) are identical). Note that "(anchor)" refers to the element of the eigenvector which has been made constant in order to make the rest of the elements of Q minimum-phase-shift. Also note that the magnitude of element (3,5) of Q is zero. Let us now summarize the main aspects discussed in this section: - The direct computation of the eigenvalues and eigenvectors of YZ with standard algorithms, such as the ones implemented in DCEIGN, is not well suited to generate the smooth functions of frequency needed to synthesize the elements of Q with rational functions. The main difficulties are eigenvector switching and closely packed eigenvalues. - These difficulties are solved by reformulating the eigenproblem from non-symmetrical to symmetrical form, and by using an extension of the Jacobi algorithm for complex symmetric matrices, in combination with good initial estimates of the eigenvector matrix (seeding) at each frequency step. - The computational speeds of DCEIGN and seeded JACOBI are essentially the same for the relatively small matrices usually encountered in cable models. - Symmetry within the elements of Q is exploited, whenever possible, in order to reduce the number of functions to be approximated. - The approximations obtained are remarkably accurate. In the examples shown, errors do not exceed 4% in the frequency ranges where the magnitudes of the functions to be approximated have practical significance. 32 1 . 2 - , F r e q u e n c y ( H z ) Fig. 7: Magnitude of the elements of column 1 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. Fig. 8: Relative error function: Magnitude, column I of Q. Evaluated using seeded JACOBI. 33 (3.1) 0 0 - 6 0 H -a " i i i n i i i n i i n i i n i I I I i inir 'L ' i i iTr i i 1 1 "! i i n n i " i n n i ' i i i i ' i r n n n , 10 10 10 1 10 10 10 10 10 10 10 Frequency (Hz) Fig. 9: Phase angle (in degrees) of the elements of column I of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. o CD (1) in (0 x; c E "o u " 5 . 0 | III l|""L I III ||""| I III II11"! I III l l " " l I III l l " " l I III l l " " l I III l l " " l I III 11 ""I I III II"'1! II,,,, , 10 10 10 1 10 10 10 10 10 10 10 Frequency (Hz) Fig. 10: Absolute error function: Phase angle (in degrees), column 1 of Q. Evaluated using seeded JACOBI. 34 CD T3 3 C an co 6" <*-o c E o 1 .2-1 1 . 0 0 . 8 0 . 6 -0 . 4 -0 . 2 -0 (4.2) (2,2) and (6.2) (3.2) (1,2) and (5.2) ^ i III i i " ^ i i l l i III n""i i III n""i i in i i ' ^ i i in Tv'<p\ Mi i i 1 1 ^ i in i i ' ^ i i in 1 1 i in n 7 10 10 10 1 10 10 10 F r e q u e n c y ( H z ) 10 10 10 10 Fig. 11: Magnitude of the elements of column 2 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. F r e q u e n c y ( H z ) Fig. 12: Phase angle (in degrees) of the elements of column 2 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. 35 1 .2-1 Frequency (Hz) Fig. 13: Magnitude of the elements of column 3 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. Frequency (Hz) Fig. 14: Phase angle (in degrees) of the elements of column 3 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. 36 1 . 2 x i 1 . 0 3 0 . 8 0 . 6 Or o 0 . 4 c 6 I 0.2H (4.4) (2.4) tnd (6.4) (3.4) (1.4) u d (S.4) " " i m i ir 1"!, i in 11""I i in i in 11""i i III i r a T T i i I in 11""i i in i in i in n , 10 10 10 1 10' 10 10 10 10 10 10 Frequency (Hz) Fig. 15: Magnitude of the elements of column 4 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. 180 GO 0) xs <v w tD X! CL Of «*-< o •<* c o u 150 120 -90 i 60 -30 0 (2.4) and (6.4) - 3 0 l""i i HI ii""L i in ii""i i HI ii""i I in 11 ""I I HI n ""i I in ii""i I in n ""i I in ii"''i I in 11 ""i i in II 7 10 10 10 1 10 10 10 10 10 10 10 Frequency (Hz) Fig. 16: Phase angle (in degrees) of the elements of column 4 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. 37 1 . 2 - 1 F r e q u e n c y ( H z ) Fig. 17: Magnitude of the elements of column 5 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. F r e q u e n c y ( H z ) Fig. 18: Phase angle (in degrees) of the elements of column 5 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. 38 .2-1 •S 1.0 o CO £ E "o 2? 0 -8 -cd 6* 0 . 6 0 . 4 0 . 2 0 (2.6) and (6.6) (1.6) and (S.6) ^ (3.6) = (4.6) = 0 X l ill I l""L l ill l ill 11 ll 11 ""l l ill l ill i l H 5 T l l l 11 ""I I III 11 ""I l ill 11 ""l i HI ii , 10 10 10 1 10 10 10 10 10 10 10 Frequency (Hz) Fig. 19: Magnitude of the elements of column 6 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. Of 180 120 -60 CO 0) 2 6 0 H si CL CD - 6 0 -c E * - 1 2 0 H o u •180 (5.6) „(1.6) (2.6) (6.6) X I ill ii""!, I ill ii""l i in ii""l I in ii""l i MI 11 ""i I ill ii""i i ill ii""i I ill 11 • • "i I ill 11 ""i i ill II , 10 10 10 1 10 10 10 10 10 10 10 Frequency (Hz) Fig. 20: Phase angle (in degrees) of the elements of column 6 of Q. Evaluated using seeded JACOBI. Exact function. Synthesized function. 39 1.4 Synthesis of the Modal Characteristic Admittance Matrix The modal characteristic admittance matrix is given by Y' = (/TP)" 1 D (36) c v yr y v 1 Y'c = ( / D ? ~ 1 Q" 1 Y Q" T- (37) Note that Y', unlike the modal propagation function A', is not uniquely defined because it depends on the normalization factors of the columns of Q. However, since these normalization factors are constrained by the minimum-phase-shift requirement imposed by the curve fitting algorithm, there is no ambiguity in the evaluation of Y ,^. In the case of underground cables, a normalization scheme that produces minimum-phase-shift functions for the elements of Q also results in minimum-phase-shift functions for the elements of Y ' . 1 2 3 4 5 6 Y' c 22 9 21 23 22 9 A' 19 14 20 12 19 16 Table 3: Order of the synthesized elements of Y' and A'. Cable length = 10 km. Figure 21 shows the magnitude of element 2 of Y' and its approximation by rational functions for the underground cable discussed above. Figure 22 shows the relative error between both curves. Figures 23 and 24 show the phase angle and error function for the same element The magnitudes and phase angles of the rest of the elements of Y' are shown in Figures 25 to 30. Note that the errors in the approximation of element 2 are larger than the errors in the approximation of the rest of the elements of 40 Y'c for this particular cable. Table 3 shows the order of the synthesized functions used. The curves shown in Figures 25 to 30 are smooth functions of frequency. Note, however, that the smoothness of the elements of Y' c depends directly on the smoothness of the elements of Q. In other words, discontinuities (or eigenvector switchovers) would reflect directly on the shape of Y ' . Fig. 2 1 : Magnitude (in S) of element 2 of Y'. Exact Junction. c Synthesized junction. g 5.0 -S 3.0 -j Frequency (Hz) Fig. 22 : Relative error junction: Magnitude of element 2 of Y'. 42 Fig. 23: Phase angle (in degrees) of element 2 of Y'. Exact function. c Synthesized junction. ^ 5.0 -i <W „ -T3 4.0 -Frequency (Hz) Fig. 24: Absolute error junction: Phase angle (in degrees) oj element 2 of Y'. CD TJ 3 •J C OO (0 O c E w F r e q u e n c y ( H z ) Fig. 25: Magnitude (in S); of elements I and 2 of Y'. Exact function. c Synthesized function. DC CD CD C/) CO JZ CL . O >-«*-< o CM c CD E CD 5 5 i in i i " ^ i i n i in 11""i i 10 10 10 1 10 10 F r e q u e n c y ( H z ) III l l " " l I III I 1 ""I 1 III II"»"I I III I I ""I I III l l"" l I I I , . . s 2 10 3 10 io5 io6 I I III II 10 Fig. 26: Phase angle (in degrees) of elements 1 and 2 of Y'. Exact function. Synthesized function. a, 0.07^ I 0.06-j CO : 5S 0.05i *~ 0.04^ 0.034 r " i i III i III i I I I i i III 11""i i III 11""i i i n 11""i i i l l 11""i i i l l 11""i i I I I i u rn 10 10 10 1 10 10 10 10 10 10 10 F r e q u e n c y ( H z ) Fig. 27: Magnitude (in S) of elements 3 and 4 of Y'. Exact function. c Synthesized Junction. F r e q u e n c y ( H z ) Fig. 28: Phase angle (in degrees) of elements 3 and 4 of Y'. Exact function. Synthesized Junction. Fig. 29: Magnitude (in S) of elements 5 and 6 of Y'. Exact function. 0 Synthesized function. CD -a w cC X! a. .o >-o CD «« lO w c a> s -5 i in 11>M>I i 10" 10 10 ""I I III II"!,'! I III l l '^ ' l I I 1 10 10 10 Frequency (Hz) M " i i in • i• • "i i in n""i 10 10 10 10' Fig. 30: Phase angle (in degrees) of elements 5 and 6 of Y'. Exact junction. Synthesized function. 46 1.5 Synthesis of the Modal Propagation Matrix Unlike the elements of Y' the elements of the modal propagation function A' are uniquely defined, regardless of the normalization factors of Q. The elements of A' are not, per se, minimum-phase-shift functions. If the magnitude of a given element of A' (e.g., A'n) is approximated with a rational function P , there is generally a good match in their magnitudes, but their phase angles will be quite different. However, as indicated in [5], there is a simple relationship between the phase angles of A' and P v e n n arg{A;i = arg{Pn}-a>rn, (39) where is defined as the time-delay constant associated with mode "n". Consequently, P n exp{-jcjT^} will match both the magnitude and the phase angle of A'. The phase angle function of the elements of A' changes very rapidly with frequency as the termor^ in equation (39) becomes dominant (see Figure 31). Therefore, rather than showing phase comparisons between A'n and its approximating function, it is preferable to show the comparison between the phase angle of P^ and the phase angle of A'n after the time-delay term has been extracted; that is, between arg{ P^ } and arglAM + CJT N. Figure 34 shows the phase angle of element 6 of A'(after COT has been extracted), and Figure 32 shows its corresponding magnitude function. Relative magnitude errors and phase differences between exact and synthesized functions are shown in Figures 33 and 35, respectively. Figures 36 to 41 show the magnitudes and phase angles for the rest of the elements of A'. Note that the phase angles corresponding to magnitudes below 0.01 have not been plotted since they are not significant from a practical point of view. Finally, Table 3 shows the order of the synthesized functions used. F r e q u e n c y ( H z ) Fig. 31: Phase angle (in degrees) of element 6 of A'. 1 .2-, T3 F r e q u e n c y ( H z ) Fig. 32: Magnitude of element 6 of A'. Exact function. Synthesized junction. 48 &T 5 . 0 -i F r e q u e n c y (Hz ) Fig. 33: Relative error junction: Magnitude of element 6 of A'. oo CD •2-<C I 0) w x: OH o CO H-> c <u £ S 30 0 - 3 0 -- 6 0 -- 9 0 •120 -•150 180 1 HIM""!, 1 III I 111 11II 1 I III I III II1"1! I III I HI 11 ""I I III ll""l I III II 7 10 10 10 1 10 10 10 10 10 10 10 F r e q u e n c y ( H z ) Fig. 34: Phase angle (in degrees) of element 6 of A' ajter extracting COT. Exact junction. Synthesized Junction. 49 5 . 0 0) -a u o u u CO CD C/l CO X! O CD .»-> C £ H - 5 . 0 -• 1 0 . 0 - 1 5 . 0 -- 2 0 . 0 "1 I III 11""! I III N""l I III 11 ""I I III II ""I I III 11 "''I I III 11 ""I I III l l " " l I III l l " " l I III l l " " l I III II 7 10 ' 10 10 1 10 10 10 10 10 10 10 Frequency (Hz) Fig. 35: Absolute error junction: Phase angle (in degrees) of element 6 of A'. The time-delay constant r determines the phase displacement between the modal propagation function and its approximation by rational functions: if the order of the approximation changes, T will also change. Also, T is numerically close to the travel time of the fastest frequency component of a wave propagating on a given mode. This can be visualized by noting that a'(t) is the response of a line or cable to an impulse 6(t); since A' has the form A'(cj) = P(a>) expl-jcjrl and a'(t) = p(t-r) u ( t - T ) , we can see that the first non-zero values of an impulse function injected at the sending end at time zero, will arrive at the receiving end at time t=r. 50 If the approximation of A' by P were exact, then r would indeed be equal to the physical travel time at very high frequency. In practice, FIT-S calculates T so that the average error in the phase angle of the approximating function is minimized. It should be mentioned that there have been earlier attempts to use FIT-S (actually an older version of the algorithm used in this thesis) to approximate the elements of A and Y c in phase quantities [16]. This approach has the obvious advantage of not requiring the approximation of Q with its inherent difficulties; that is, the solution of Y V +1 = F = A F, cm m m k Y V, - I, = B, = A B c k k k m would not require the transition between modal and phase quantities. However, early in this thesis project, it was found that the elements of A (phase) could not be approximated directly with acceptable accuracy. In part, this is due to the fact that the magnitude of the elements of A can be highly oscillatory, as shown in Figures 42 and 43. Also, even if the magnitude of the elements of A could be matched successfully, their phase angles cannot be matched using a single time-delay term. In [16], the approximating functions were passed right through the maxima and minima of the oscillatory parts of A, and the discrepancies in the phase angles were either ignored or simply dismissed as unimportant It could be argued that it is unnecessary to fit the magnitudes of the elements of Y', A' and Q at frequencies as high as 10 MHz, when the maximum frequency which can be reproduced in a numerical solution is limited by the sampling rate At (f /71ajc= l/2At). It would not be wise, however, to tie the computationally intensive procedure of fitting the cable transfer functions to the time step of a particular transient 51 simulation, given that the same cable (i.e., its transfer functions in synthesized form) can be used in many situations, each of which may require a different value of At Also, from a database-management point of view, it is desirable to keep the generation of cable parameters and the transient simulations as separate procedures. Recent work carried out by J. R. Marti at UBC has resulted in a "post-processing" algorithm to reduce the order of the approximations of A' and Y' according to the time step used in the transient simulation. By "post-processing", we mean that the order reduction is carried out in the EMTP, and not in the cable constants support program. This permits the storage of a high-order, wide-range fit in the database, and the use of a reduced-order fit in the time simulation. Computational savings of up to 30% in the time-step loop of the EMTP can be obtained with this procedure (with no loss of accuracy). Let us now summarize the main points covered in this section: - In order to perform the numerical convolutions needed to update the past history current sources of the equivalent circuit of the new cable model, the elements of Y', A' and Q must be synthesized with rational functions in the frequency domain. To this effect, a curve-fitting algorithm (FIT-S) has been used. To obtain a good match in both magnitude and phase angle, the functions to be approximated must be smooth, minimum-phase-shift functions of frequency. - Standard eigenvalue/eigenvector algorithms do not produce smoothly varying eigenvectors as a function of frequency. To avoid this difficulty, a new eigenvalue/eigenvector algorithm, based on Jacobi's method for real symmetric matrices, has been developed for the complex symmetric case. Also, instead of solving the unsymmetric eigenproblem 52 YZ Q = Q the following equivalent problem is solved (assuming Y and Z in loop quantities): H R = R D , yz' where H = / T Z v/T and the modal transformation matrix Q is given by Q = / T R . The eigenvectors obtained with this method are smooth functions of frequency which can be synthesized very accurately with FIT-S. - Once the difficulties in the evaluation of Q are solved, fitting the elements of Y' c and A' is relatively straightforward: If the elements of Q are smooth functions of frequency, then the elements of Y' will also be smooth functions of frequency. Unlike the elements of A in phase quantities, the elements of A' can be synthesized with rational functions because it is possible to associate a single time-delay term with each propagation mode. V T3 3 C CO vs b-i o w CD 1 . 2 '"(. I III ll""l I III ll""l I III I III I III ll""l I III 11 ""I I III ll""l I III ll""l 10 10 10 1 10 10 10 10 10 F r e q u e n c y ( H z ) 10 10 Fig. 36: Magnitude of elements 1 and 2 of A'. Exact Junction. Synthesized Junction. tuO CD T3 CO E -•S I <L> w CO 0-O CV? •a 7 1 H - > c S 30 0 - 3 0 4 - 6 0 - 9 0 -3 - 1 2 0 - 1 5 0 -E - 1 8 0 -i - 2 1 0 - 2 4 0 ^ - 2 7 0 \ i in ii""i I in ii""i i in 11""i I in ii""i I in 10 10 10 1 10 10' ,1 i ill 11""i i III ii""i i III ii""i i ill n""i i in II 7 2 10 3 1 0 4 10 5 10 6 10 7 F r e q u e n c y ( H z ) Fig. 37: Phase angle (in degrees) of elements 1 and 2 of A'ajler extracting COT. Exact Junction. Synthesized Junction. 54 "j. 111111 " " i I in I in i i ' 10 10 10 1 10 10 10 Frequency (Hz) Fig. 38: Magnitude of elements 3 and 4 of A'. Exact function. Synthesized function. oo CD T3 3 cd H )e I CD m cO o CO OT c CD e CD 3 30 0 -30 -. -60 -i -90 -j -120 -j -150 -= -180 -210 -240 -270 4111111-^ 1 10 10 rrr™ 10" 1 imii""i nun""! iinn"''i iiuii""i iiini""i iiini""i n u n , 10 10 10 10 10 10 10 Frequency (Hz) Fig. 39: Phase angle (in degrees) of elements 3 and 4 of A' after extracting UT. Exact function. Synthesized function. CD X i 3 1 . 2-1 "c 1 . 0 OC 2 o co in w c CL) E CD 3 0 . 8 -0 . 6 -0 . 4 0 . 2 0 "1 I III I I III 11 ""I I III 11 ""I I III 11 ""I I III 11 "''I ! Ill 11 ""I I 10 10 10 1 10 10 10 F r e q u e n c y ( H z ) Fig. 40: Magnitude of elements 5 and 6 of A'. Exact function. Synthesized function. 11inn"j.'i i i l ' i i T ^ i i n m 10 10 10 oo <D X! 3 CO •s I <D CO co CL 3 30 -q 0 - 3 0 -. - 6 0 -j - 9 0 \ -< "120 - j o - 1 5 0 \ co -- 1 8 0 -E <*S --«o - 2 1 0 - j CO : E " 2 7 0 ^- ^ i in n " ^ i in 111:1! i in i i " " i i in 11""i i in \ \™\ i in i i ' ^ ' i i in ii"j;i i III i i ' ^ i i III \\™\ i III IK, 10 1 0 ' 10 10 10 10 1 . 1 0 10 10' F r e q u e n c y ( H z ) Fig. 41: Phase angle (in degrees) of elements 5 and 6 of A' after extracting ur Exact function. Synthesized function. 10 56 1 . 2 - j CD Frequency (Hz) Fig. 43: Magnitude of element (2,2) of A (phase). 57 1.6 Diagonalization of YZ In the derivation of the equations which led to the equivalent circuit of the new cable model, two major assumptions were made: i) The elements of q(t), y'(t) and a'(t) can be expressed accurately as finite sums of C exponentials. ii) The matrix product YZ is diagonalizable; that is, it is possible to find a matrix Q such that Q" 1 YZ Q = D . yz The first assumption was justified in the previous sections, where it was shown that the elements of Q, Y^ and A' can be accurately synthesized with rational functions using FIT-S. We will now address the second assumption; that is, whether it is always possible to find the eigenvector matrix Q. First, it should be noted that the unsymmetric eigenproblem YZ Q = Q D is equivalent to the modified, symmetric problem H R = R D . yi where Q = / T R as long as Y and Z are expressed in loop quantities. If Y and Z are expressed in node quantities, then equations (28) to (31) should be used instead. 58 If Q is invertible, then R will also be invertible because in the case of underground cables, Y (in loop quantities) is diagonal with nonzero elements; therefore, if Q" 1 exists then R" 1 will also exist, and viceversa. We emphasize this point because the theoretical considerations made in this section will refer to the symmetric problem in order to take advantage of the special properties inherent to symmetric matrices. From a mathematical point of view, the existence of a set of linearly independent eigenvectors of a symmetric matrix H is assured as long as all the eigenvalues of H are distinct Difficulties arise when the multiplicity of one or more of the eigenvalues is greater than one (e.g., two or more identical eigenvalues). In this case, a set of linearly independent eigenvectors may or may not exist [18]. This situation created a certain amount of controversy for a number of years, where mathematicians claimed that a modal transformation matrix could not always be found, while engineers unfailingly found modal transformation matrices in their practical applications [19]. This controversy was recently resolved (at least for the case of overhead transmission lines) by J. A. Brandao Faria and J. F. Borges da Silva [20], who isolated a set of physical conditions (tower configuration, earth resistivity, etc.) where a modal transformation matrix could not be found at a specific frequency in the presence of two identical eigenvalues. From a practical point of view, however, the situation where two or more identical eigenvalues may cause a failure in the evaluation of the eigenvector matrix can be avoided completely without significant loss of accuracy. A. J. Gruodis and C. S. Chang [21] proposed disturbing the elements of the diagonal elements of YZ with random noise of very small magnitude (about four orders of magnitude smaller than the magnitude of the diagonal entries). This approach effectively reduces to zero the probability of finding numerically identical eigenvalues, thus guaranteeing the existence of 59 Q at all frequencies. This disturbance does not have a significant effect on the numerical values of the eigenvalues. In this thesis, it was found that when the modified eigenproblem is solved using seeded JACOBI, eigenvector instability due to closely packed (or even identical eigenvalues) does not occur, even in cases where DCEIGN produces erratic results. This can be explained by the way the modified Jacobi algorithm finds the eigenvalues and eigenvectors of H. If H is seeded with the eigenvector matrix obtained in the previous frequency step (i.e., H, = R5 1 H R0), the matrix to solve, H , , is practically diagonal already. Since with JACOBI, the eigenvectors of are found by small rotations of R, (which is practically the unit matrix), numerical instability does not occur. Additionally, the seeding procedure introduces a small, random numerical disturbance in because of the additional arithmetic needed. In summary, - From a mathematical point of view, it is possible to encounter situations where multiple eigenvalues do not permit the evaluation of the modal transformation matrix Q at a given frequency. - In practice, however, the probability of encountering an irregular eigenvector problem can be reduced to zero when the diagonal elements of the matrix to be diagonalized are disturbed randomly by a small amount - A similar effect is obtained when the eigenproblem is solved in modified form using seeded JACOBI. Consequently, from a numerical point of view, the assumption that YZ is diagonalizable at all frequencies can be accepted as valid. 60 CHAPTER 2 Numerical Results In order to establish the accuracy of the new cable model, four types of simulations will be used: 1) Open and short-circuit responses in the frequency domain. 2) Transient simulations where the new cable model is compared with exact analytical solutions. 3) Transient simulations comparing the response of the new cable model with the response of the EMTP's frequency dependence line model with constant transformation matrix. 4) Comparison with a field test With the first two types of simulation, it is possible to determine, theoretically, how accurate the new cable model is when the input data (in this case the cable parameters) are assumed to be correct The comparison with the state-of-the-art line model of the EMTP, illustrates the range of validity and the limitations of existing time-domain line models. The comparison with the results of the field test adds a degree of confidence to the purely theoretical results obtained with the theoretical simulations. The first three theoretical comparisons are performed on the 230 kV underground cable described in Appendix m. Physical data for the 110 kV crossbonded cable used in the field test comparison is also given in Appendix III. 61 2.1 Open and Short-Circuit Responses in the Frequency Domain Consider the single-phase transmission system shown in Figure 44. ///////////////////////////////////////////////////// Fig. 44: Single-phase transmission system. If the receiving end is open, the receiving-end voltage will be given by m ' 2 A'(w) Vo 1 + A' 2(CJ) (40) On the other hand, if the receiving end is short-circuited, the receiving end current \ ' will be given by 2 Y ; M A ' ( O > ) V 0 I'(w) = -M 1 - A' 2(GJ) (41) Equations (40) and (41) represent the open and short-circuit responses of a transmission system, at a given frequency u> [22]. Note that the sending-end voltage source V© is constant Primed symbols have been used to indicate that the single-phase system shown in Figure 44 can also represent a single propagation mode of a multi-conductor system. The extension to phase quantities is then made through the relationships 62 (42) = Q - T M V '(CJ). (43) where A'(w) and Y^(CJ) now represent diagonal matrices, and Vo is a vector which depends on the actual voltage sources (in phase quantities); that is, Therefore, open and short-circuit responses of an underground cable can be obtained (at any given frequency) from equations (40) to (43) when Y^ , A' and Q are known. The exact frequency response is obtained when Y', A' and Q are used directly. The response of the new model is obtained when these matrices are calculated from their respective approximations by rational functions. To obtain the response of the frequency dependence line model with constant transformation matrix5 currently used in the EMTP, matrices Y'c, A' and Q are calculated as follows: 1) The transformation matrix Q is evaluated at a suitable frequency CJ 0 • The columns of Q(CJ 0) are then rotated so that the imaginary parts of their elements are minimized in the least-squares sense [23]. 2) The imaginary part of Q(CJ 0) is set to zero; the resulting matrix Q 0 = Re{Q(w0)}, is then assumed to be the transformation matrix which diagonalizes Y Z over the entire frequency range. Therefore, Y' and A' are calculated (at all frequencies) using Qo instead of Q(o>). 3) The elements of Y' and A' are synthesized with rational functions. Finally, Q 0 and the approximations by rational functions of Y' and A' are used to calculate the open 'Future references to the "frequency dependence line model" should be understood as references to the "frequency dependence line model with constant transformation matrix". Vo = V 0 Q T (c4 63 and short-circuit responses of the frequency dependence line model. The results obtained in the frequency domain with equations (40) to (43) can also be duplicated in the time domain using the EMTP. A transient simulation, which starts from steady-state initial conditions (at co=co0), and where the system is not disturbed after time zero, produces a sinusoidal response of frequency co0. Note that the peak value of this response corresponds to the real part of the frequency-domain response at co = co0. Comparisons between these time-domain simulations and the direct frequency-domain calculations show that both answers are virtually identical, as long as the time step of the transient simulation is suitably small (e.g., At ^ 0.2n7 co0). Therefore, open and short-circuit responses in the frequency domain provide an indirect way to measure the accuracy of a model under transient conditions; that is, a good agreement in the frequency domain indicates that the transient response will also be good. In the following frequency-domain simulations, two different forms of excitation will be used: a) Symmetric excitation b) Zero-sequence excitation. 64 a) Symmetric Excitation For this set of open and short-circuit responses, the sheaths of the reference three-phase cable are grounded at the sending end (see Figure 45); the magnitudes of the voltage sources connected to the cores are set to 1.0 V, and their phase angles are set 120° apart (i.e., -120°, 0° and +120°). Note that this form of excitation is analogous to the case where a three-phase transmission line is connected to balanced sources with three ground wires grounded at the sending end. ///////////////////////////////////////////////////// Fig. 45: Three-phase underground cable. Symmetric excitation. Figure 46 shows the magnitude of the open-circuit response (i.e., the receiving-end voltage for core 1) of the new cable model. Figures 47, 48 and 49 show the response of the frequency dependence line model, when Q is evaluated at 0.1 Hz, 60 Hz and 5 kHz, respectively. Each plot contains two superimposed curves, with the exact response shown in dashed trace. Only the receiving-end voltage for core 1 is shown because it is representative of the response of the other conductors. Figure 50 shows the magnitude of the short-circuit response (i.e., the receiving-end current for core 1) of the new model. Figures 51, 52 and 53 show the short-circuit response of the frequency dependence line model, where Q has also been 65 evaluated at 0.1 Hz, 60 Hz and 5 kHz, respectively. From these results, it can be seen that the open and short-circuit responses of the new cable model are accurate over the entire frequency range considered. Also, the open-circuit response of the frequency dependence line model with Q evaluated at 5 kHz is remarkably good. However, the short-circuit response with Q evaluated at 5 kHz is only accurate at frequencies above 1 kHz. At frequencies below this point, the errors are substantial. Note, that when Q is evaluated at 60 Hz, the short-circuit response does not match the exact response at this frequency. The reason for this discrepancy is that the imaginary part of Q at 60 Hz (which is set to zero when Y' and A' are calculated) is too large to be neglected. Also, note that the phase angles of the open and short-circuit responses follow the same behaviour of their magnitudes; that is, if the magnitudes are close to the exact analytical response, then the phase angles will also be close (see Figures 54 and 55). 66 Fig. 46: Open-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. New model. cO cr CD GO CO o > •a c cu I to c CD o 0) 9 . 0 8 . 0 -7 . 0 : 6 . 0 ~ 5 . 0 -j 4 . 0 -f 3 . 0 2 . 0 ^ 1 . 0 11II 11""i 11II 11""i 11II 11""i 111 10" 1 10 10 11 " , " l I III I l ' " ' l I I I I 11''"I I III 111 'J'l T U M I 10 10 10 10 10 F r e q u e n c y ( H z ) Fig. 47: Open-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. Constant Q, evaluated at 0.1 Hz. N X O CD CO cy CU Ofl cO > c CD I c CD CJ CD 7 . 0 n 6 . 0 5 . 0 4 . 0 3 . 0 -j 2 . 0 -j 1 . 0 - • "I l l l l l l " " ! I l l l l l " " ! M l l l l " , " l I M i l l " " ! 1 III 11'' *T1 M I M I ' ' " I I I I M T^'fTl II 10"" 1 10 10 ' 1 0 ' 10" F r e q u e n c y ( H z ) 10 10 10 Fig. 48: Open-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. Constant Q, evaluated at 60 Hz. S 7 . O n F r e q u e n c y ( H z ) £ 4 . 0 n 1 0 " F r e q u e n c y ( H z ) Fig. 49: Open-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. Constant Q, evaluated at 5 kHz. Fig. 50: Short-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. New model. Fig. 5 1 : Short-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. Constant Q, evaluated at O.J Hz. 0 I'' " i — i i l l I i 1 ' " i — i i l l i I 1 1 1 11—I i l l i I " " i — i i l l i i " " i — i i l l II 10" 2 4 10" 2 4 10" 2 4 1 2 4 10 2 4 10 Frequency (Hz) Fig. 52: Short-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. Constant Q, evaluated at 60 Hz. Fig. 53: Short-circuit response (magnitude). Symmetric excitation: Core 1. Exact response. Constant Q, evaluated at 5 kHz. 73 o c CD t-u. 3 CJ XI c d 0) oCD Frequency (Hz) Fig. 54: Short-circuit response (phase angle). Symmetric excitation: Core 1. Exact response. New model. m or c CD u u 3 o x i C CD BO c 0) o m OS • 2 7 0 | | | | [ | 1 1 1 '| | | | | | | M i i | | | || | | i n i | I I | | i | i n i | | | || | | | | | | | » | | H I I , 10 1 10 10 10 10 10 10 Frequency (Hz) Fig. 55: Short-circuit response (phase angle). Symmetric excitation: Core 1. Exact response. Constant Q, evaluated at 5 kHz. 74 b) Zero-Sequence Excitation For the following set of simulations, a single voltage source of amplitude 1.0 V and zero phase angle is connected to all six conductors at the sending end of the cable (see Figure 56). This form of excitation is analogous to the excitation of the zero-sequence mode of a balanced six-phase transmission line. -o--o c, -O S , —o c2 -o s2 -o-—o c3 O S 3 ///////////////////////////////////////////////////// Fig. 56: Three-phase underground cable. Zero-sequence excitation. Figure 57 shows the open-circuit response (receiving-end voltage for core 1) of the new model. Figures 58, 59 and 60 show the response of the frequency dependence line model when Q is evaluated at 0.1 Hz, 60 Hz and 5 kHz, respectively. Figures 61 to 64 show the corresponding short-circuit responses (receiving-end currents for core 1). The same general comments made in the case of balanced excitation are also valid here. Note, however, that the short-circuit response when Q is evaluated at 5 kHz is considerably poorer than its counterpart for the balanced case. In fact, the receiving-end current for core 1 is practically zero at all frequencies (see Figure 64). In summary, - The open and short-circuit responses of the new cable model are very accurate over the entire frequency range of interest 75 - The frequency dependence line model with constant transformation matrix produces acceptable answers over a limited frequency range, as long as the frequency at which Q is evaluated is properly chosen; that is, as long as the imaginary part of its elements are small. This is generally the case when Q is evaluated at very high, or at very low frequencies. U 5.0 -1 O s F r e q u e n c y ( H z ) F r e q u e n c y ( H z ) Fig. 57: Open-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. New model. 5 . 0 - i 4 . 0 -3.0 -2 . 0 -. 0 ""I I III 11 • • • • I M l l l l " ' 1 ! I III I r ' j T T T T n f 'L"l I M i l l 1 ' " ! M I N I 10 ' 1 10 10 10 10 10 10 10 Frequency (Hz) 58: Open-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. Constant Q, evaluated at O.I Hz. £ 5 . 0 . o CD Io 4 . 0 -Frequency (Hz) 59: Open-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. Constant Q, evaluated at 60 Hz. 78 Frequency (Hz) Fig. 60: Open-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. Constant Q, evaluated at 5 kHz. Fig. 61: Short-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. New model. Fig. 62: Short-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. Constant Q, evaluated at 0.1 Hz. 1 2 . 0 1 0 . 0 -8 . 0 6 . 0 H 4 . 0 2 . 0 0 l 1 ' " i i III i i 1 ' " i i i i i i i — i i i i i i " " i II 10 2 4 10 2 4 10" 2 4 1 2 4 Frequency (Hz) rm-\ i 10 2 4 10 0 . 10 0 . 0 9 0 . 0 8 0 . 0 7 0 . 0 6 0 . 0 5 0 . 0 4 0 . 0 3 0 . 0 2 0 . 0 1 0 1 ' 1 ' | " " | " " | " l " l " l ' V n ' 1 1 i |Mn|'"i|M|"|'VVV1 ' i ' i p m p i t p T p i p ^ - r n M 1 1 " l " " l ' 'I"I"I'VI^ 10 2 3 4 6 10 2 3 4 6 10 2 3 4 6 10 2 3 4 6 10 Frequency (Hz) Fig. 63: Short-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. Constant Q, evaluated at 60 Hz. 82 0 . 1 0 0 . 0 9 0 . 0 8 0 . 0 7 0 . 0 6 0 . 0 5 0 . 0 4 0 . 0 3 0 . 0 2 0 . 0 1 0 10 T p m p u 111' <' |""I' '11 'I'Tl'VI i , n u | ! i n | i i r 1 i l , i V V ^ , , , , | i M , | U M | i r l|>V1'Vl , , , , , 5 2 3 4 6 I f / 2 3 4 6 10 2 3 4 6 10 2 3 4 6 10 F r e q u e n c y (Hz ) ac 1 2 . 0 1 0 . 0 8 . 0 5 6 . 0 4 . 0 -ao 2 . 0 -c T T T T l — I I I I I I ' ' " I—I I I I l I " . " l— I I I I I I " "I—I I I I I I " " I — l l ITT1 , 10 2 4 10 2 4 10 2 4 1 2 4 10 2 4 10 F r e q u e n c y ( H z ) Fig. 64: Short-circuit response (magnitude). Zero-sequence excitation: Core 1. Exact response. Constant Q, evaluated at 5 kHz. 83 2.2 Impulse and Square-Wave Responses in the Time Domain Consider the three-phase underground cable shown in Figure 65, where a voltage source v (t) is connected to core 1. Cores 2 and 3, as well as all the sheaths, are grounded at the sending end. v, (0 6 — o c, -o s, -o c 2 -o s, — o c 3 -o s 3 ///////////////////////////////////////////////////// Fig. 65: Three-phase, 230 kV underground cable. If v0(t) is a sinusoidal source of frequency w 0, the exact steady-state response of this cable system can be found analytically (e.g., by solving equations (3) and (4) in phasor form at O> = C J 0 ) . If v0(t) is not sinusoidal, it can be approximated by a Fourier series of sine and cosine terms; that is, vo(0 * ygenes = ao+ 2 & n cos(nco0t) + b sin(nu0t)], n= I (44) where OJ0 = 2ir/T, and T is the time period over which v0(t) is approximated by ^series'6 e x a c t r e s P o n s e to this series can be evaluated by superimposing the steady-state responses of every term of the expansion; that is, the response at frequencies CL>=ncj0, where n=l,2,...,N, and N tends to infinity. 'If v0(t) is periodic, then the Fourier series becomes an exact representation of v0(t) in the limit, as n tends to infinity. If v0(t) is not periodic, then the Fourier series expansion only approximates v0(t) for 0 ^ t ^ T. 84 From a numerical point of view, it is not possible to express v with an S6ri€S infinite number of terms; therefore, the exact analytical response to an arbitrary input waveform v0(t) cannot be be evaluated without a certain amount of error. However, if a sufficiently large number of terms "N" is taken into account, the resulting finite series, N vo(t) - v (t) = a0+ £ [a cos(ncj0t) + b n sin(ncj0t)], (45) <J . - m i l fl 1=1 will be a reasonable approximation of v0(t) for O^t^T. Furthermore, the exact analytical response of v^ (t) can be readily obtained using superposition. From these considerations, we now have the means to obtain an analytically exact response (in the time domain) of an n-phase cable which can be used as a reference solution to assess the accuracy of the transient response of the new cable model. a) Impulse Response Figure 66 shows the finite Fourier series approximation v^ (t), of a 1.2 x S.Ojxs voltage impulse v0(t) of the form v0(t) = k [expC-t/T,}- exp{-t/T2l] (46) with T, = 3.48MS, T 2 = 0.80MS, and k = 2.014 V. This series approximation contains 500 terms, and the period T is 1 ms. The situation simulated is the energization of the reference 230 kV, three-phase underground cable, with v^ (t) as shown in Figure 65 (the length of the cable is assumed to be 1 km). The exact response is calculated using superposition; that is, by adding up the steady-state response of every term of the Fourier series approximation of v0(t). The response of the the new cable model. is obtained using the EMTP. Note that in the 85 EMTP simulation, v^ (t) (rather than v0(t)) is used as input, in order to be consistent with the analytical calculations. Figure 67 shows the voltage at the receiving end of. core 1. The analytical solution is indicated with a dashed line, while the EMTP simulation is shown in solid trace. Note that the Fourier series approximation of v0(t) is a periodic function; therefore, the response will also be periodic Since the transient has not died completely after one period T, the analytical solution does not start from zero initial conditions when the impulse is injected at times 0, 1 ms, 2 ms, etc. The EMTP solution, on the other hand, does start from zero initial conditions at t=0. This is the reason for the discrepancies in the first cycle. In subsequent cycles, the EMTP simulation converges to the exact solution as proper initial conditions are met Figure 68 shows the third cycle of the voltage at the receiving end of core 1, and Figure 69 shows the voltage at the receiving end of sheath 1 for the same time frame. Note that the response of the new cable model and the exact analytical solution are virtually identical at this point of the simulation. 1 . 10 T 1 . 0 0 0 . 9 0 -3 0 . 8 0 0 . 7 0 0 . 6 0 ^ 0 . 5 0 0 . 4 0 0 . 3 0 -3 0 . 2 0 0.10-3 0 - 0 . 1 0 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 .0 1 . 2 1 . 4 T i m e ( m i l l i s e c o n d s ) 1 . 6 1 . 8 2 . 0 2 . 2 Fig. 66: Fourier series approximation of a voltage impulse. N=500, T=l ms. 87 - | — i — i — I — | — I — i — I — | — I — i — I — [ — i — i — I — ) — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — [ — ' ' r -0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 1 . 2 1 . 4 1 . 6 1 . 8 2 . 0 2 . 2 T i m e ( m i l l i s e c o n d s ) Fig. 67: Impulse response, receiving-end voltage^ Core 1: Analytical response. New model. 88 1 . 0 0 - i - 1 — I — j — i — I — i — i — | — t — r — i — i — | — I — I — I — I — | — r — r ~ I — I — [ — i — I — I — I — | — i — i — i — i — | — I — I — i — i — | 2 . 0 2 . 1 2 . 2 2 . 3 2 . 4 2 . 5 2 . 6 2 . 7 2 . 8 2 . 9 3 . 0 T i m e ( m i l l i s e c o n d s ) Fig. 68: Impulse response, receiving-end voltage (third cycle). Core 1: Analytical response. New model. — 0 . 4 0 ^ CD * 0 . 3 0 H 0 . 2 0 -5 o.io CO CD s: 0 8 -0.10-3 c o w - 0 . 2 0 CD I-CD - 0 . 3 0 CO ct - 0 . 4 0 1—i—i—i—p—i—i—i—r—|—i—r-T "! "]—i—i—t r - j T -1—i—i—pi i i i — | — i i t i j l — i — I i | i—i—i—t j i i — i — i | £ 2 . 0 2 . 1 2 . 2 2 . 3 2 . 4 2 . 5 2 . 6 2 . 7 2 . 8 2 . 9 3 . 0 T i m e ( m i l l i s e c o n d s ) Fig. 69: Impulse response, receiving-end voltage (third cycle). Sheath 1: Analytical response. New model. 89 b) Square-Wave Response In this simulation, the voltage source in Figure 65 now represents a square-wave generator. The Fourier series approximation of the input voltage contains 1000 terms; and the period T is 40 ms (see Figure 70). The length of the cable is assumed to be 10 km. Figure 71 shows the receiving-end voltage for core 1, and Figure 73 shows the receiving-end voltage for sheath 1. Note that the discrepancies between analytical and numerical solutions are considerably smaller than in the case of the impulse response. Since the analytical solution starts from nearly-zero initial conditions (at times 0, 40 ms, 80 ms, etc.), the initial error in the EMTP simulation is very small. After the artificial transient caused by the discrepancy in the initial conditions has attenuated, analytical and numerical solutions become practically' identical (see Figure 72). The transient response of the cable used in this case is considerably attenuated after the first 20 ms of the simulation. This is the reason why the analytical solution starts from nearly-zero initial conditions. Therefore, during the first half of each cycle, the square-wave response also represents a reasonable approximation of the step response of this cable. 1 .60 1 .40 i 1 .20 1 .00 1 0 .80 -j 0.60 i 0 .40 0 .20 -31 0 -0.20 4 -0 .40 -0.60 0 -T—i—i—r- I—I—I—T 10 20 30 40 50 60 T i m e ( m i l l i s e c o n d s ) 70 80 90 -0 . .10 -q .00 \ .90 \ .80 \ .70 -j .60 4 .50 | .40 -j .30 \ .20 \ : 10 \ 0 -j-10 L i i i I I i r - r " | l T 1—r—|—i i t i j i— r n — I j i -i r ~ i—p - i — i T T y r — 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 T i m e ( m i l l i s e c o n d s ) —i—|—i—i—i—r 4.0 4.5 5. Fig. 70: Fourier series approximation of a square wave. N=1000, T=40 ms. CD GO CO -4-> > CD !_ O O CD c O D-, to CD CD > CO I CD U CO 3 cr 1 . 6 0 1 . 4 0 -3 1 . 2 0 1 . 0 0 0 . 8 0 0 . 6 0 0 . 4 0 0 . 2 0 -3 0 - 0 . 2 0 - 0 . 4 0 i - 0 . 6 0 i l I l 0 10 20 30 4 0 5 0 60 T i m e ( m i l l i s e c o n d s ) 70 8 0 90 CD GO CO -t-> o > CD t, O O CD CO C o a w CD U CD > CO £ I CD U o 3 a* w 1 . 6 0 -q 1 . 4 0 ^ 1 . 2 0 1 . 0 0 0 . 8 0 -3 0 . 6 0 0 . 4 0 0 . 2 0 -3 0 - 0 . 2 0 0 ~ l— l— I — I — I — I — I — I — r - - l — I — I — | — l — I " ~ 1 — 1 — I — I — I— I— I — T " 6 8 10 12 14 T i m e ( m i l l i s e c o n d s ) 16 18 2 0 Fig. 71: Square-wave response, receiving-end voltage. Core 1: Analytical response. New model. CO oo CO "o > CO s-o o CO V) c o a OT CO u CO > tO is I CO u co 3 cr CO .60 n .40 .20 -j .00 -j .80 .60 -3 .40 .20 - 0 . 2 0 40 42 -j—i—(—i—i—|—i—i—T~T—j—i—i—i—(—j—i 44 46 48 50 1 ' I ' ' ' ' l ' 46 48 50 52 54 T i m e ( m i l l i s e c o n d s ) -i—»—i—i—n—i—r~ 56 58 60 Fig. 72: Square-wave response, receiving-end voltage (second cycle). Core 1: Analytical response. New model. M CO "o > CO CO sz ai c o a VI CO t. > CO I CO CO 3 cr CO -i—i—i—i—i—i—i—i—r—]—i—i—i—r—j—i—i—i—i—|—i—i—i—r—|—r—i—i—i—|—i—i—i—i | i—i—i—i—|—i—i—i—i—| 42 44 46 48 50 52 54 56 58 60 T i m e ( m i l l i s e c o n d s ) Fig. 73: Square-wave response, receiving-end voltage (second cycle). Sheath I: Analytical response. New model. 93 2.3 Simulation of a Single-Phase Line-to-Ground Fault The effects of taking into account the variation with frequency of the modal transformation matrix Q were highlighted in section 2.1 by means of open and short-circuit responses in the frequency domain. These simulations gave a quantitative measurement of the accuracy of the steady-state response of a model over the entire frequency range of interest From these results we know, for example, that the short-circuit response of the frequency dependence line model with Q evaluated at 5 kHz has an error of the order of 30% in the vicinity of 200 Hz, under balanced excitation. However, this information cannot be related directly to the behaviour of the model during a transient simulation. In other words, we know that the accuracy of the response of the frequency dependence line model is limited to a certain frequency range, but it is not evident how different its transient response would be when compared to the transient response of the new cable model. To illustrate the effect of taking into account the frequency dependence of Q during a transient simulation, a single-phase line-to-ground fault will be simulated with the EMTP. The connection diagram for this simulation is shown in Figure 74. —o 6 c, —^ o o s, <^ >\—o o c2 — ' o o s2 ///////////////////////////////////////////////////// t = 10 ms Fig. 74: Simulation of a single-phase line-to-ground fault. 94 The peak magnitude of the 60 Hz voltage sources is 1.0 V, and their phase angles are set 120° apart (i.e., -120°, 0° and +120°). The cores of the 230 kV reference cable are open at the receiving end; and the length of the cable is assumed to be 10 km. The simulation starts from 60 Hz steady-state initial conditions; and the receiving-end of core 3 is connected to ground at t=10ms. Figures 75, 76 and 77 show the receiving-end voltages when the sheaths are grounded at both ends of the cable. Figure 78 shows the current at the receiving end of the faulted phase. The response of the new cable model is shown in solid trace; and the response of the frequency dependence line model with Q evaluated at 5 kHz is shown in dashed trace. Figures 79 to 82 show the results of the same simulation when the sheaths are open at both ends. Note that the receiving-end voltages of both models are identical during the first 10 ms of the simulation. This is to be expected, given the open-circuit responses of these models in the frequency domain (see Figures 46 and 49). When the sheaths are grounded, the overvoltages at the receiving end of the unfaulted phases are considerably reduced because the circulating currents in the sheaths oppose any sudden voltage changes. Therefore, the differences between the fault currents shown in Figure 78, are very similar to the differences observed in the short-circuit responses at 60 Hz with balanced excitation (see Figures 50 and 53). When the sheaths are not grounded, there is no transient-suppression effect Therefore, it is practically impossible to establish a direct link between the frequency-domain responses and the transient simulation. Nevertheless, these results illustrate the importance of taking into account the frequency dependence of Q in the simulation of transients in underground cables. Fig. 75: Receiving-end voltage. Core 3 (sheaths are grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. 2.0-1 £ 1 .5-o u o —2 .0 —1—•—'—>—i—'—'—'—'—r-1—1—1—1—i—1—1—1—1—i—1—1—1—1—i—1—1—1—1—i—1—1—1—1—r~1—1—'—1—i 0 5 10 15 20 25 30 35 40 Time (mi l l i seconds) Fig. 76: Receiving-end voltage. Core 2 (sheaths are grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. Fig. 77: Receiving-end voltage. Core 1 (sheaths are grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. 0.5 -. —0 .6 —'— 1— i—<—j— 1 —'—i—i— i—'—•— i—'— i—'—»—'— i—j—'— 1— 1—'— i— 1 —'—>—'—i—'—*—'— 1— j—•— 1—•—•— i 0 5 10 15 20 25 30 35 40 Time (milliseconds) Fig. 78: Receiving-end current. Core 3 (sheaths are grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. 97 2.0 n CO c '.H - 1 . 5 -o : OJ —2.0 —1—•—'—<—j—'—•—>—'—i—'—<—'—*—j—i—*—'—'—j—'—'—'—>—|—i—i—i—'—|—i—i—i—i—]—i—<—'—'—) 0 5 10 15 20 25 30 35 40 Time (milliseconds) Fig. 79: Receiving-end voltage. Core 3 (sheaths are not grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. —2 .0 —'—1—•—1—i—•—'—1—'—i—•—'—1—'—i—>—•—•—'—i—>—1—'—'—i—•—•—>—•—i—1—1—•—r~~\—•—'—•—'—i 0 5 10 15 20 25 30 35 40 Time (milliseconds) Fig. 80: Receiving-end voltage. Core 2 (sheaths are not grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. Fig. 81: Receiving-end voltage. Core J (sheaths are not grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. 0 .20 -a co o c u 3 o 3 (0 10 15 20 25 Time (milliseconds) 30 35 Fig. 82: Receiving-end current. Core 3 (sheaths are not grounded): Frequency dependence line model with Q evaluated at 5 kHz. New model. 99 2.4 Comparison with a Field Test The field test presented here has been reproduced from [9], which in turn quotes [24] as the original source. Figure 83 shows the circuit diagram of the test, and Appendix Ul summarizes the physical data of the crossbonded cable used. 500 a soon Fig. 83: Field test connection diagram. A voltage impulse of waveshape Ox40MS (i.e., with negligible virtual front time) and a peak magnitude of 7.3 kV was applied between the centre core and earth. The sheaths are connected together and grounded through a 10 0 resistance at the sending and receiving ends. Figure 84 shows the measured and the calculated voltages at the sending end of core 2. Figures 85 to 87 show the voltages at the second crossbonding point Measured results are shown in dashed trace; and the calculated results are shown in solid trace. From these graphs it can be seen that measured and calculated results agree reasonably well. Given that [9] questions the accuracy of the cable data, a second simulation was performed In this simulation, the virtual time to half-value of the input voltage was 100 changed from 40 ns. to 45 M S . The sending-end (core 2) voltage for the second simulation is shown in Figure 88. Examination of this curve shows an improved agreement with the experimental results. However, the calculated sheath voltages at the second crossbonding point were not affected by this change. In a third simulation, the earth resistivity was reduced from 100 Om to 50 Om. In this case the agreement between calculated and experimental results also improved to a certain extent The discrepancies between measured and calculated results could be attributed to several factors: a) Insufficient data regarding dielectric losses and their variation with frequency. b) Sheath voltages were measured with respect to ground. Sheath to sheath measurements would have been more reliable, and the effect of equipment grounding would not have been significant c) The description of the input waveform may have been oversimplified or loosely described (for example, a virtual front time of 0 /is is not physically possible), as suggested by the improvement in the calculated results when the virtual time to half-value was changed by 10%. Overall, the agreement between the waveshapes of calculated and simulated results are reasonably good, considering the many uncertainties involved in such a test 101 1 .50 1 1 1 i 1 1 1 i ' 1 1 i 1 ' • i • 1 1 i 1 1 <-y 0 20 40 60 80 100 120 140 160 180 200 T i m e ( m i c r o s e c o n d s ) Fig. 84: Sending-end voltage. Core 2 (input waveshape 0x40us): field test. New model. 0.15 I""" I — I — | — i — i — i — r ™ ~ i — I — i — | — r — T — l — | — i — r — ? — | — i " I — I — ] — l — I — I — j — i — i — 1 — | — i — i — i — r — r — r — i — | 0 20 40 60 80 100 120 140 160 180 200 T i m e ( m i c r o s e c o n d s ) Fig. 85: Second crossbonding-point voltage. Sheath I (input waveshape 0x40ns): field test. New model. 102 0.15-1 ~i—|—i—i—i—|—i—i—i—|—i—i—i—|—i—i—i—|—i—i—i—|—i—i—i—|—i—r—1—| 0 20 40 60 80 100 120 140 160 180 200 Time (microseconds) Fig. 86: Second crossbonding-point voltage. Sheath 2 (input waveshape 0x40n%): field test. New model. 0 20 40 60 80 100 120 140 160 180 200 Time (microseconds) Fig. 87: Second crossbonding-point voltage. Sheath 3 (input waveshape 0x40n%): field test. New model. 103 1 .50 - i :> --•—• • 1 .20 GO <C _ .*_) -> 0 -90 -C\J -CD _ u o _ o 0 6 0 -1J -X) _ o 2 0 3 0 --<D -• - I — I — I — I — I — I — I — I — I — I — 1 — I — I — I — I — I — 1 — I — I — I — I — I — I — [ -0 20 40 60 80 100 120 140 160 Time (microseconds) i ' 1 1 i 180 200 Fig. 88: Sending-end voltage. Core 2 (input waveshape 0x45us): field test. New model. 104 2.5 Numerical Stability From an analytical point of view, the stability of the model is assured by the fact that all the poles in the closed-form approximations of Y' A'and Q are negative; therefore, the exponential terms exp { -pi} cannot run away during a time-domain simulation. However, for overall numerical stability, the stability of the integration rule used in the discretization of the differential equations is required as well. In the derivation of the equivalent circuit for the new model, the type of integration rule used appears implicitly in the coefficients of the recursive convolutions indicated in Appendix II. These coefficients are known as the "Recursive Convolution" integration coefficients, and they produce completely stable solutions [25]. In the solution algorithm of the EMTP, on the other hand, the "Trapezoidal" rule of integration is used; and this integration rule is also stable.7 Consequently, the implementation of the new cable model into a host program such as the EMTP produces stable numerical solutions. An extra measure of stability is provided by the fact that Q" 1 is not approximated directly by rational functions. From equation (19) it can be seen that, to obtain the modal currents, the following numerical convolution must be evaluated i'(t) = q- '(t) • i (t). However, i'(t) is obtained indirectly; that is, m i (t) = q(t) • i'(t) i (0 = qoi'(0 + h' (t) i'(t) = q61 i (t) - q51 h' (t). 7 Under some circumstances, the Trapezoidal rule of integration produces numerical oscillations. However, the amplitude of these oscillations is always bounded [26]. 105 From a numerical point of view, this means that Q Q" 1 = U will be true regardless of the errors incurred during the fitting process. In turn, this means that the new cable model is "internally consistent"; that is, errors in the synthesized forms of Y', A' and Q, would only be apparent if the original Y and Z matrices were to be reconstructed. If we call these reconstructed matrices Y 0 and Z 0 , the relationships Q - 1 exp{-i/To"Zo7} Q = A' (47) Q A'Q" 1 = exp{-/YoZo7} (48) will always be true, and free of fitting errors. If, on the other hand, Q" 1 were to be approximated as well, the inevitable errors in the fitting process would prevent equations (47) and (48) from being strictly true. This would result in an accumulation of errors during a transient simulation. Note that if both Q and Q" 1 were approximated with rational functions, the resulting model would also be stable, but the accumulation of errors might restrict its accuracy to relatively short transient simulations. In practice, the above reasoning was supported by the fact that the simulations shown in section 2.2 were allowed to proceed for relatively long simulation times (in excess of 12000 time steps) and that no deviations from the analytical solutions were detected. 106 2.6 Computational Speed In this thesis project, the main concern has been to obtain a very accurate' cable model capable of reproducing the frequency dependence of Y^, A' and Q. In other words, emphasis was centered on accuracy and not on computational speed, since the new model was to take into account what no other EMTP line or cable model could. As long as the computational speed had not been several orders of magnitude lower than that of existing models, the new model would have been considered useful from a practical point of view. It was, then, a most rewarding surprise to discover that the computational speed of the new model was very high. To measure the relative computational speed of the new cable model, the energization of the reference 230 kV cable was simulated using UBC's version of the EMTP. No network components, other than the voltage sources, were included in the simulation of this six-conductor system. Measurements of computational speed were based on the elapsed CPU time during the time-step loop of the EMTP. Two of the most commonly used EMTP line models were considered in this comparison: i) Constant-parameter line model, where the line is assumed to be balanced and lossless; the series resistance of the line is taken into account by lumping it in three places. ii) Frequency dependence line model, where the frequency dependence of Y' and A' is taken into account, but the transformation matrix Q is assumed to be constant. For this test, Q was chosen to be the well-known a,/3,o transformation matrix. The order of the approximations by rational functions of the elements of Y' c and A' was forced to be the same as the ones used in the new cable model. 107 Model Relative timing Constant- parameter 1.00 Frequency dependence 6.00 New model. 6.14 Table 4: Relative CPU times. EMTP time-step loop. Table 4 shows the results of these comparisons. Note that the speed of the new cable model is practically the same as the speed of the frequency dependence model with constant Q. This is remarkable, given the additional computational effort required to perform matrix-vector convolutions rather than matrix-vector products. Although these measurements are probably influenced by the way the different models have been programmed, they indicate that the additional computational overhead needed to take into account the frequency dependence of Q is not prohibitively high. 108 CHAPTER 3 Future Research in Cable Modelling and Related Areas The research effort involved in this thesis project has answered a number of questions regarding the modelling of underground cables in the time domain, especially regarding the behaviour of their modal transformation matrices. However, as more knowledge is gained on a subject, new questions arise, and the need for further research becomes apparent What follows is a summary of the main topics which, in the author's opinion, will require additional research in the future. 3.1 Simulation of Pipe-Type Underground Cables The number of underground cables studied in this thesis project was, by necessity, relatively small. Although underground cables of similar construction should behave along the same general lines, other types of cables, such as pipe-type cables, may present different characteristics in the behaviour of Q. Interfacing existing pipe-type cable software with the cable constants support program developed in this thesis project would probably be the easiest way to examine the behaviour of the transformation matrices for this type of cables. 3.2 Extension of the New Cable Model to Overhead Transmission Lines The cable model presented in this thesis should be, in theory, directly applicable to any n-conductor system. Therefore, it could be used in the simulation of multiple-circuit overhead transmission lines. Double-circuit lines, for instance, are 109 notoriously difficult to simulate accurately with existing line models which assume constant Q, when a broad range of frequencies is involved in the transient simulation1 Some work on the subject has already been done as part of this thesis project; however, it is not reported here because time limitations did not permit the presentation of well-documented results. The main difficulty in the application of the new cable model to multiple-circuit transmission lines lies in the behaviour of the modal transformation matrix Q. The magnitude of the elements of Q associated with most of the double-circuit lines studied at UBC changed smoothly with frequency (see Figure 89). However, it appears that in cases of strong asymmetry (e.g., when a 155 kV circuit and a 230 kV circuit share the same tower) the magnitude of the elements of Q may change very sharply with frequency, as shown in Figure 90. Analysis of the eigenvectors of one of these asymmetrical cases indicated that these variations are not caused by the eigenvector switchover phenomena discussed in Chapter 2. These peaks are not a major problem because FIT-S can be modified to introduce complex conjugate poles capable of following such magnitude functions. However, unlike the eigenvectors of underground cables, it appears that the eigenvectors of overhead transmission lines cannot be reduced to minimum-phase-shift functions by means of simple scaling. Although the accurate reproduction of the magnitude of the elements of Q would represent a considerable improvement over assuming Q to be constant, it is clear that the accuracy achieved in the case of underground cables cannot be duplicated in the case of transmission lines, as long as the phase angles of the elements of Q are not reproduced accurately. •For example, in the simulation of fast transients, when a good match for 60 Hz steady-state initial conditions is also required. 110 Future research on this problem could be focused on several alternatives. For example, it has been found that successive approximation, in additive form, improves the approximation of the phase functions of the elements of Q for overhead lines; that is, if a given function H0is first approximated with H , , and H 0 - is approximated with H 2 , the agreement between arg { H, + H, } and arg { H 0 1 is much better than the agreement between argtH, land arg{ H 0 1 . Another possibility is the use of all-pass functions to compensate for the differences between arglHjland arg{H0l. Ill 1 . 0 0 01 7J 0 . 8 0 -4-1 C GO <0 S 0 . 6 0 tir c 6 3 0 . 2 0 o o ••"i IIIIM""J l l l l l l " " l 1 111 II' • • »1 l l l l l l " ' ' l l l l l l l " " l I III 11 ""I I I I IH""I I III II 10" 10" 1 10 10 10 10 10 10 10 F r e q u e n c y ( H z ) 89: Magnitude of column 1 of Q of a double-circuit transmission line. Two 500 kV circuits sharing the same tower. 2 . 0 -1 . 8 — CD _ ud 1 6 --i-> -c 1 .A _ oo -a 1 2 — Of 1 0 -««-< _ o 0 8 _ c 0 6 — 6 -3 0 A --o 0 2 0 -10 10 10 '""I I III H""l I III 11 ""I I III l l " l ' l I III 11 ""I I III l l " " l I III l l " " l I III II 7 ^ io3 io io 5 io6 io 7 1 10 10 F r e q u e n c y ( H z ) 90: Magnitude of column 1 of Q of a double-circuit transmission line. 115 kV and 230 kV circuits sharing the same tower. 112 CONCLUSIONS A new model to simulate electromagnetic transients in underground high-voltage cables has been presented in this thesis. This new model belongs to the class of time-domain models, and it is directly compatible with the solution algorithm of the EMTP. However, unlike existing time-domain models, the frequency dependence of the modal transformation matrix Q is accurately taken into account Conceptually, the new cable model is simple: the variation with frequency of Y', A' and Q is represented in closed form by means of rational functions. This implies that the behaviour of these transfer functions in the time domain can be expressed by finite sums of exponential functions; this representation, in turn, permits the recursive solution of the differential equations which describe the transient behaviour of an underground cable. One of the main difficulties encountered was the evaluation of the elements of the modal transformation Q as smooth functions of frequency, so that they could be approximated with rational functions. Since standard eigenvalue/eigenvector algorithms are ill-suited for this purpose, a new strategy for the solution of the eigenvalues and eigenvectors of the matrix product YZ was developed. The new algorithm first converts the eigenproblem from unsymmetric to symmetric form. The eigenvalues and eigenvectors are then calculated with an extension of the Jacobi method (modified to handle complex symmetric matrices). Continuity between the eigenvectors calculated at two adjacent frequency points is obtained by seeding; that is, by using the eigenvector matrix obtained in the previous frequency step as an initial guess in the evaluation of the eigenvector matrix. The eigenvalues and eigenvectors thus obtained are the smooth functions of frequency needed for the synthesis of the elements of Y ' , A' and Q by rational 113 functions. The accuracy of the synthesized forms of Y', A" and Q is very high. This accuracy is generally higher than the accuracy with which earth resistivity, dielectric losses, and other cable parameters can be estimated. Therefore, from a practical point of view, the new cable model is only limited by the accuracy of the input data, and by the simulation capabilities of its host program. Computational speed is very high. Therefore, the detailed simulation of crossbonded cables using the new model becomes attractive, from a practical point of view. Also, given that minor sections are usually of the same length, only the parameters of a single section need to be synthesized. Furthermore, by modelling each section separately (making the crossbonding connections explicitly), non-linear voltage limiters at the crossbonding points can be easily taken into account The new model is numerically stable. In relatively long simulations (of the order of 12000 time steps) there has been no noticeable accumulation of errors. The new model is internally consistent; that is, the transition between phase and modal quantities (and viceversa), is not affected by the accuracy of the approximation by rational functions of the elements of Q. This is possible because the equations leading to the equivalent circuit of the new cable model have been manipulated so that the approximation by rational functions of the elements of Q" 1 is not needed explicitly. Finally, the new model presented in this thesis is general; its extension to the simulation of multiple-circuit overhead transmission lines could also be of considerable practical importance. 114 APPENDIX I Evaluation of Cable Parameters A program to calculate the cable parameters required as input data for the new cable model has been written as part of this thesis project This program calculates the parameters for an underground cable system consisting of n single-phase cables. Each cable has two or more metallic conductors, of which, one is the central "core" conductor, and the others are concentric sheaths or armours (see Figure 1.1). The formulae used are based on the following assumptions: a) The cable system consists of conductors whose axes are mutually parallel, and are parallel to the surface of the earth. b) The system is longitudinally homogeneous. In particular, the earth is assumed to be homogeneous, laterally and longitudinally. c) Displacement currents are assumed to be negligible. d) Proximity and end-effects are ignored. A brief summary of the formulae used in the cable constants support program is given below. LA Shunt Admittance and Series Impedance Matrices in Loop Quantities Consider a cable system consisting of n single-phase cables. Each single-phase unit consists of m coaxial conductors (e.g., core, sheath and armour, for m=3), as shown in Figure 1.1. 115 earth insulation insulation insulation core sheath armour Fig. 1.1: Cross- section of an underground cable. The shunt admittance matrix can be formulated more easily if loop equations are used; that is, d I dx loop _ y V loop loop 0.1) where loop Y l l 0 0 0 0 0 Y[2 0 (1.2) ll3 *loop m^ ^loop 3 1 6 ^ e v0^2^t and current vectors (of dimension nm) in loop quantities, (see Figure L2); Y^, Y ^ and Y ^ are submatrices dimensioned mxm. Y^, for instance, has the form 116 f //////////////////////////////// nl Vi ///////////////////////////////////. Fig. L2: Single-phase underground cable. 11 cs 0 0 0 0 0 Y 0 (1.3) where Y is the admittance between core and sheath, cs Y is the admittance between sheath and armour, sa Y a e is the admittance between armour and ground, and Y ^ for example, is given by Y = G + jcoC . cs cs J cf (1.4) G = G +o>C tan6, cs ocs cs (1.5) where Ccs is the capacitance in F/km, between core and sheath; G£s is the 117 corresponding shunt conductance (G is the shunt conductance at dc); and 6 is defined ocs as the loss angle (typical values of tan6 for modern cables range around 0.002). Formulae for the evaluation of these capacitances can be found in [27]. Even though the loss angle 8 depends on frequency [28], in the cable constants support program, tan5 is assumed to be constant because data on the variation of the loss angle with frequency is very difficult to obtain. Note that if all single-phase units are identical, then = Y ^ Similarly, for the series impedance matrix Z ^ ^ , we have d V, ^ = Z. I. (1.6) d x loop loop ' "loop Z Z Z in nn ni3 Z121 Z122 Z123 Z131 Z132 Z133 (1.7) where Z^ is symmetric and dimensioned n m x n m; Z^. is a symmetric submatrix dimensioned mxm. A diagonal submatrix of Z^ ^ o r Stance) has the form 777 Z Z 0 c cs Z Z z cs s sa z z sa a (1-8) where 118 Z £ is the impedance of the loop consisting of the core with return through the sheath. Zs is the impedance of the loop consisting of the sheath with return through the armour. Z f l is the impedance of the loop consisting of the armour with return through the earth. Zcs is the mutual impedance between the core-sheath and sheath-armour loops. Z ^ is the mutual impedance between the sheath-armour and earth return loops. Note that there is no coupling between the core-sheath and armour-earth loops. The off-diagonal submatrices of Z^ have the form 0 0 0 0 0 0 0 0 eij (1.9) where Z .. is the mutual earth return impedance between cables i and j. Note that only BIJ the armour-earth loops of each single-phase unit are coupled. Also note, that if the cables are identical, the diagonal submatrices in (1.7) will also be identical. The formulae used to evaluate the impedances which do not include earth return terms were obtained from [27], while the formulae for the earth return impedances were obtained from [29]. LB Transformation Between Loop and Phase Quantities From Figure 1.2 it can be seen that for a single-phase cable, the transformation between loop and phase voltages and currents can be derived from 119 V. = V . - V . 11 nl n2 Y12 = Yn2 " Yn3 Y13 = Yn3 nl 11 lnl ~ ll2 ~ lll ln3 = ll3 " ll2 • (1.10) From (1.10), matrices Z and Y in phase quantities can be obtained through manipulation of the rows and columns of Z ^ ^ and Y^ [27]. However, a more general way to make the transformation, would be to use the relationships V, = A V , loop phase (I.H) I, = B I , , loop phase (1.12) where A T = B" 1 , and U - u 0 u 0 0 A = 0 u - u B = u u 0 (1.13) 0 0 u u u u where U in equation (1.13) represents an identity submatrix of dimension (mxm). Introducing (1.11) and (1.12) into (1.1) and (1.6), results in the following relationship: Y. Z, = BY , Z . B" 1 loop loop phase phase (1.14) Note that B and A are orthonormal matrices, and that pre and post-multiplication °^ ^phaseZphase ^ ^ ^ * n v e r s e ' represents a similarity transformation. Therefore, YloopZloop ^ ^phaseZphase ^ a v e ^ e s a m e e*8envalues. The relationship between the eigenvectors in loop and phase quantities is given by 120 Q, = BQ , . (1.15) loop phase ' In the cable constants support program, the evaluation of the eigenvalues and eigenvectors is made in loop quantities. Since Y'c and A'are both the same in loop and phase quantities, only Q must be transformed into phase form before its elements are approximated by rational functions. I.C The Generalized Eigenproblem Consider the unsymmetric eigenproblem Y Z Q = QD, (1.16) where D is a diagonal matrix containing the eigenvalues of YZ, and Q is the corresponding eigenvector matrix. Let S be the eigenvector matrix which diagonalizes the shunt admittance matrix Y; then, S- 1 Y S = Dy (1.17) where is a diagonal matrix containing the eigenvalues of Y. Introducing (1.17) into (1.16) gives SD S ' 1 Z Q = QD y •D^S- 1 ZQ = (/D^)- 1 S" 1 QD /D^S- 1 Z [S (/TT)(/D^- 1 S" M Q = (/D^" 1 S" 1 Q D (1.18) 121 Let us define H and R as H = / T X S - 1 Z S / T T (1.19) R = (/D^)- 1 S" 1 Q. (1.20) Introducing these definitions into (1.18) gives HR = RD (1.21) S / T T R = Q. (1.22) Equation (1.21) represents the eigenproblem associated with matrix H. Note the following: a) Matrix S (which diagonalizes Y) satisfies the condition S" 1 = ST [17], and matrix H in equation (1.19) is symmetric. b) Since H is symmetric, the eigenvector matrix R also satisfies R" 1 = R T . c) The eigenvalues of H are the same eigenvalues of YZ. d) If the eigenvectors of H are known, the eigenvectors of YZ can be obtained using equation (1.22). As indicated in Appendix LA, Yj is a diagonal matrix; therefore, S = U, and the modified eigenproblem is reduced to H. R. = R. D (1.23) loop loop loop v ' ^loop loop ^ loop^^loop (1-24) loop ^ loop ^loop (1-25) 122 Some of the advantages of formulating the eigenproblem in generalized form using loop equations are: a) The eigenvector matrix R is orthonormal; therefore, no matrix inversions are needed to evaluate either R"1 or Q" 1 (see equation (1.25)). b) H is symmetric; therefore, the modified Jacobi method can be used to find D and R. c) With the use of appropiate seeding, this procedure is faster than conventional methods (for the type of cables studied in this thesis). It should be noted that all mathematical considerations concerning the eigenvector matrices Q, R and S, depend on the assumption that YZ, H and Y, respectively, are diagonalizable. As it was indicated in section 1.6, these matrices are diagonalizable as long as their eigenvalues are distinct Even when eigenvalues coalesce, it is possible to obtain eigenvector matrices using the perturbation techniques suggested in [21]. However, from a mathematical point of view, the eigenvectors resulting from quasi-identical eigenvalues may or may not be ill-conditioned. Therefore, the validity of the numerical solutions should be constrained by the existence of uniformly bounded condition numbers [31]; that is, KQ) = IQU ||Q" 1EU * L (independent of u). 123 APPENDIX H Numerical Recursive Convolutions Let q(t) be a finite sum of exponentials of the form m p(t) = k0 6(t-T) + I k exp{-p/t-r)} for t£ T , and q(t) = 0 for t£ T. Then the numerical convolution of q(t) with a given function f(t) can be expressed as g(t) = q(t) • f(t) m g(t) = k0f(t-r) + I g/t), (II.1) /=/ where, g/t) = c.g/t-r) + d^t - T ) + e^t-r-At), and c^ , d^ , and e^ . are integration constants which depend on k f p^ . and At [5]. If q(t) is the time-domain representation of the approximation by rational functions of the elements of Y' or Q, then T = 0, and (II.1) becomes g(t) = df(t) + h(t), where m d = ko+ I d/t) 1=1 1 m h(t) = I c. g/t-At) + e f(t-At). 124 Note that h(t) depends only on past history values of f and g. If q(t) corresponds to the approximation of the elements of A', then k0 = 0, and (II.l) becomes m g(t) = I c .^g/t-At) + d.f(t-r) + e.f(t-T-At). /=/ Therefore, g(t) depends exclusively on past history values of g^ and f. 125 APPENDIX m Physical Data of Test Cables Physical data for the 230 kV underground cable used in all the theoretical comparisons of this thesis (Cable 1), and for the 110 kV crossbonded cable used in the comparison with the field test (Cable 2), are given below. Cable 1 Cable 2 r, (cm) 0.00 1.05 r2 (cm) 2.34 1.78 r3 (cm) 3.85 3.08 rfl (cm) 4.13 3.26 rs (cm) 4.84 4.08 Core Resistivity (flm) 0.0170 10' 6 0.0183 10-Sheath Resistivity (flm) 0.2100 10-6 0.0280 10-Inner insulation tanh 0.001 0.040 Outer insulation tanh 0.001 0.100 Inner insulation e 3.5 3.3 Outer insulation e r 8.0 3.8 Earth resistivity (Qm) 50 100 /////////////////////////////////////// •if-0 0 0 1.2 m y .25 m .25 m—y-Fig. in.l: Cable 1 configuration. /////////////////////////// 1.0 m 0 0 .35 m 0 ¥ 35 m Fig. HI.2: Cable 2 configuration. 127 REFERENCES [I] A. Semlyen and A. Dabuleanu, "Fast and Accurate Switching Transient Calculations on Transmission Lines with Ground Return using Recursive Convolutions." TKFF. Transactions on Power Apparatus and Systems, pp. 561-571, March/April 1975. [2] W. S. Meyer and H. W. Dommel, "Numerical Modelling of Frequency-Dependent Transmission-Line Parameters in an Electromagnetic Transients Program". IEEE Transactions on Power Apparatus and Systems, vol. PAS-93, pp. 1401-1409, September/October 1974. [3] J. F. Hauer, "State-Space Modelling of Transmission Line Dynamics via Nonlinear Optimization". IEEE Transactions on Power Apparatus and Systems, pp. 4918-4925, December 1981. [4] A. Ametani, "A Highly Efficient Method for Calculating Transmission Line Transients". IEEE Transactions on Power Apparatus and Systems, vol. PAS-95, No. 5, pp. 1545-1551, September/October 1976. [5] J. R. Marti, "Accurate Modelling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Calculations." IEEE Transactions on Power Apparatus and Systems, pp. 147-157, January 1982. [6] L. M. Wedepohl and C. S. Indulkar, "Switching Overvoltages in Long Crossbonded Cable Systems Using the Fourier Transform." IEEE Transactions on Power Apparatus and Systems, pp. 1476-1480, July/August 1979. [7] H. Nakanishi and A. Ametani, "Transient Calculation of a Transmission Line Using Superposition Law". IEE Proceedings, vol. 133, PL C, No 5, pp. 263-269, July 1986. [8] CIGRE Working Group 13.05, "The Calculation of Switching Surges: II.- Network Representation for Energization and Re-Energization Studies on Line Fed by an Inductive Source". Electro, No. 32, pp. 17-42, 1974. [9] N. Nagaoka and A. Ametani, "Transient Calculations on Crossbonded Cables." IEEE Transactions on Power Apparatus and Systems, April 1983, pp. 779-787. [10] H. W. Dommel and W. S. Meyer, "Computations of Electromagnetic Transients". IEEE Proceedings, vol. 62(7), pp. 983-993, July 1974. [II] A. Papoulis, The Fourier Integral and its Applications. McGraw-Hill Book Company, New York, San Francisco, Toronto, pp. 204-208, 1962. [12] B. T. Smith, et al., "Matrix Eigensystem Routines - HSPACK Guide". Lecture Notes in Computer Science, Springer, Verlag, Berlin, Heidelberg, New York, 1974. [13] Personal communication with H. W. Dommel. [14] B. N. Parlett, The Symmetric Eigenvalue Problem. Prentice-Hall Inc., pp. 309-310, 1980. [15] R. W. Hornbeck, Numerical Methods. Quantum Publishers INC, pp. 263-240, 1985.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Simulation of electromagnetic transients in underground...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Simulation of electromagnetic transients in underground cables with frequency-dependent modal transformation… Marti, Luis 1986
pdf
Page Metadata
Item Metadata
Title | Simulation of electromagnetic transients in underground cables with frequency-dependent modal transformation matrics |
Creator |
Marti, Luis |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | This thesis presents a new model to simulate the behaviour of underground cable systems under transient conditions. The new cable model belongs to the class of time-domain, frequency-dependent models, and it is directly compatible with the solution algorithm of the EMTP (Electromagnetic Transients Program). The most important feature of the new model is that it takes into account the frequency dependence of the modal transformation matrices and cable parameters, thus overcoming the main limitation of currently-used transmission line and cable models, which assume that the modal transformation matrices are constant Conceptually, the new model is relatively simple. The system parameters which define the behaviour of an underground cable (namely the modal characteristic admittance matrix, the modal propagation matrix, and the modal transformation matrix), are expressed in closed form by approximating them with rational functions in the frequency domain. Therefore, in the time domain, all numerical convolutions can be expressed recursively. The host transients program (to which the model is interfaced) sees the new model as a constant, real admittance matrix, in parallel with a continuously-updated vector current source. The accurate approximation by rational functions of the modal transformation matrix is possible when its elements are continuous and smooth functions of frequency. Standard eigenvalue/eigenvector algorithms are not well suited for this purpose. Therefore, a new procedure to generate eigenvalues and eigenvectors has been developed. This procedure is based on the Jacobi method, and it produces the desired smooth functions of frequency. This manuscript presents a number of simulations where the performance of the new cable model is compared with exact analytical solutions. These simulations show an excellent agreement between analytical and numerical answers. The effects of not taking into account the frequency dependence of the modal transformation matrices is illustrated with the simulation of a line-to-ground fault on a three-phase cable. The response of the new cable model is also compared with results measured in a field test The new cable model is numerically stable. Its computational speed is comparable to that of frequency-dependent line models with constant transformation matrices. The new cable model is general. Its extension to the simulation of multiple-circuit overhead transmission lines should also be of considerable practical importance. |
Subject |
Underground electric lines |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0064979 |
URI | http://hdl.handle.net/2429/27442 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1987_A1 M34_3.pdf [ 4.87MB ]
- Metadata
- JSON: 831-1.0064979.json
- JSON-LD: 831-1.0064979-ld.json
- RDF/XML (Pretty): 831-1.0064979-rdf.xml
- RDF/JSON: 831-1.0064979-rdf.json
- Turtle: 831-1.0064979-turtle.txt
- N-Triples: 831-1.0064979-rdf-ntriples.txt
- Original Record: 831-1.0064979-source.json
- Full Text
- 831-1.0064979-fulltext.txt
- Citation
- 831-1.0064979.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064979/manifest