C A L C U L A T I O N OF F R E Q U E N C Y - D E P E N D E N T PARAMETERS OF U N D E R G R O U N D POWER CABLES WITH FINITE E L E M E N T M E T H O D By Y A N A N YIN B. Eng., Shanghai Jiaotong University, 1982 M . Sc., Tsinghua University, 1984 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September 1990 © Yanan Yin, 1990 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 copying agree that permission for extensive of this thesis for scholarly purposes may be granted department or by his or her representatives. It is by the head of my understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia Vancouver, Canada Date DE-6 (2/88) StifMtvtU^ 26j m o Abstract In this thesis, the finite element method (FEM) is applied to the calculation of frequencydependent series impedances and shunt capacitances of underground power cables. The principal equations describing the quasi-magnetic fields and static electric fields are solved with FEM based on the Galerkin technique. The Js method and the loss-energy method are derived to calculate the impedances of a multiconductor system from its field solution, and the energy method and the surface charge method are derived to calculate the capacitances. With a single-core (SC) coaxial cable, the suitability of quadratic isoparametric elements and high-order simplex elements are studied, and a suitable division scheme is suggested for the auto-mesh program. The conventional FEM with a field truncation boundary is applied to the impedance calculation of buried SC cables. Suitable locations for the field truncation boundary and division schemes in the earth are studied. The results show that r^, > 12S is required e to obtain accurate impedances of shallowly buried cables with the conventional FEM. This requires a large solution region in the earth at low frequencies. A new technique based on the perturbation concept is proposed to reduce the solution region in the earth. Comparisons between the results from the conventional FEM and from the proposed technique with a significantly reduced solution region in the earth show good agreement. In the case studies, the FEM is applied to the parameter calculation of multiphase SC cables, PT cables, sector-shaped cables, and stranded conductors. The numerical results are compared with those from analytical formulas. ii Table of Contents Abstract ii List of Tables ix List of Figures xii Acknowledgements xiii 1 Introduction 1 2 Impedance Calculation with Finite Element Method 8 2.1 Introduction 8 2.2 Definition of Transmission Line Parameters 8 2.3 Principal Equations in Impedance Calculations 11 2.4 FEM Based on the Galerkin Technique for the Principal Equations . . . 14 2.4.1 The Galerkin technique for solving differential equations 14 2.4.2 The FEM based on the Galerkin technique 16 2.5 2.6 2.7 [Z] Calculations from the Field Solutions 20 2.5.1 The J method 20 2.5.2 The loss-energy method 21 s Impedances of Single-Core Coaxial Cables Summary 24 26 3 Element Types and Shape Functions iii 28 3.1 Introduction 28 3.2 High-Order Simplex Elements 29 3.2.1 Shape functions in simplex elements 31 3.2.2 Integral matrices of simplex elements 32 3.3 36 3.3.1 Quadratic quadrilateral isoparametric element 36 3.3.2 Quadratic triangular isoparametric element 39 3.4 Calculation of Integrals in the Loss-Energy Method 3.5 General Procedures for [Z] calculations with FEM 3.6 3.7 4 Isoparametric Elements [Z] Calculation of an SC Coaxial Cable with FEM 42 42 45 3.6.1 Optimum divisions for SC coaxial cables 47 3.6.2 Computation efficiency: CPU time, storage, and pivoting . . . . 55 Summary 58 Earth Region Reduction Technique for [Z] Calculation with F E M 4.1 Introduction 4.2 Analytical Formulas for [Z] of Buried SC Coaxial Cables 4.3 [Z] Calculation of Deeply Buried SC Coaxial Cables by Conventional 60 FEM with a Field Truncation Boundary 4.4 Earth Reduction Technique 4.5 [Z] Calculations of Shallowly Buried SC Coaxial Cables with FEM . . . 4.5.1 62 67 74 82 Determination of the solution region for the conventional FEM in [Z] calculations of shallowly buried cables 4.5.2 60 83 Application of the proposed technique to [Z] calculations of shallowly buried SC cables 86 iv 4.5.3 Comparisons between analytical results and FEM results for shallowly buried SC cables 4.5.4 4.6 91 [Z] calculations for a cable layout of arbitrary structure Summary 95 5 Admittance Calculation with Finite Element Method 97 5.1 Introduction 97 5.2 Principal Equation and FEM solution 98 5.3 5.4 5.5 6 93 [C] Calculation from the Field Solutions 100 5.3.1 The energy method 101 5.3.2 The surface charge method 103 [C] Calculation of SC Coaxial Cables 107 5.4.1 General form of [C] for SC coaxial cables 107 5.4.2 [C] calculation of a SC coaxial cable by FEM 108 Summary 112 Case Studies in [Z] and [C] Calculations 114 6.1 Introduction 6.2 [Z] Calculations of Buried or Tunnel Installed Multiphase Cable Systems 115 6.3 [Z] and [C] Calculations of PT Cables 6.4 114 119 6.3.1 The [Z] calculation of PT cables 119 6.3.2 The [C] calculation of PT cables 125 [Z] and [C] Calculations of Sector-Shaped Cables 128 6.4.1 The [Z] calculation of sector-shaped cables 128 6.4.2 The [G] calculation of sector-shaped cables 130 6.5 The Calculation of Internal Resistance of Stranded Conductors 131 6.6 Summary 135 v 7 Conclusions and Recommendations for Future Work References 137 141 A Integral matrices [Q^] and [T ] of simplex elements 146 B Detailed Derivation of Pollaczek's Formula 150 s C List of Symbols 156 vi List of Tables 3.1 [QW] index strings 35 3.2 Locations of sampling points and weighting factors for Gaussian quadrature 39 3.3 Locations of sampling points and weighting factors for quadratic triangular element 41 3.4 Division factors f&. for the SC coaxial cable 49 3.5 Division radii for iso and sim2 49 3.6 [R] and [L] of a two-conductor SC coaxial cable 50 3.7 Storage and other parameters for different elements 56 3.8 CPU time requirements for different elements 56 3.9 Storage and CPU time requirements for pivoting 57 3.10 Storage and other parameters for different elements 57 3.11 CPU time requirements for different elements 58 3.12 Storage and CPU time requirements for pivoting 58 4.1 Earth division radius patterns for isoparametric elements in the decade from 10 to 10 n 70 n+1 4.2 Earth division radii for isoparametric elements 70 4.3 [R] and [L] of the deeply buried SC coaxial cable 71 4.4 [Z] of the deeply buried cable found from a reduced earth region 4.5 [Z] found with the proposed technique based on Ep 81 4.6 [Z] of the shallowly buried SC coaxial cable from the conventional FEM . 84 4.7 [Z] of the shallowly buried SC coaxial cable from the proposed technique 87 vii . . . . 79 4.8 I . . . . 87 4.9 Storage and other parameters for the proposed technique (rfc=5m) . . . . 89 ep found with numerical integrations for the proposed technique 4.10 CPU time requirements for the proposed technique (rt=5m) 89 4.11 Storage and other parameters for the conventional FEM (rj,=12£ ) . . . . 89 4.12 CPU time requirements for the conventional FEM (rj,=12£ ) 90 e e 4.13 Maximum differences in [Z] with the proposed technique at different p e . 90 e . 90 . . . 91 4.16 Threshold / with Pollaczek's formula for maximum differences <1% . . . 92 4.17 Threshold / with Pollaczek's formula for maximum differences< 10% . . 92 4.18 Threshold / with Pollaczek's formula for maximum differences<30% 92 4.14 Maximum differences in [Z] with the proposed technique at different r 4.15 Maximum differences in [Z] with Pollaczek's formula at different p e 4.19 [Z] of the tunnel installed cable from three approaches . . 94 5.1 i ando inP (A ,0= iEfc o C 5.2 Radius divisions in [C] calculation of the SC coaxial cable 109 5.3 [C] of a two-conductor SC coaxial cable . 110 5.4 Surface |E| for different divisions with isoparametric elements Ill 6.1 [^AAK[^CCD 6.2 [•ZABK^BCD an< 6.3 [^AAKI^CCD a n 6.4 [•^ABK^BC]) 6.5 [^AAI 6.6 [•ZABK^ACD 6.7 [^AAKI^CCI) 6.8 [Z^g]([Zg(j]) and [Z^Q] of the PT cable with a cradle arrangement . . 126 6.9 [C] of the PT cable with the triangle arrangement (/iF/km) r m m i A N ( m p A N < ^ [^BBI 0 °^ d [^BB] °^ * n e [^ACl °^ ^ [^BBK^CC}) an< [^BBI °^ 117 buried three-phase cable system 117 tunnel installed three-phase cable system 118 * °^ ^ [^BCl °f * ied three-phase cable system Dur ^ [^ACl °^ u n n e ^ installed three-phase cable system 118 c a n e 106 < m i ca ^ with a triangle arrangement . 123 e ble with a triangle arrangement . 123 PT cable with a cradle arrangement . . 126 viii 127 6.10 [C] of the PT cable with the cradle arrangement (/zF/km) 127 6.11 [Z] of the sector-shaped cable 129 6.12 [C] of the sector-shaped cable (/zF/km) 131 6.13 The internal resistance of the two-layer stranded conductor 133 6.14 The internal resistance of the one-layer stranded conductor 134 6.15 The internal resistance of the three-layer stranded conductor 135 6.16 The internal resistance of the four-layer stranded conductor 135 A.l [Q^] and [Ts] of the 1st order simplex element A.2 [Q^\ and [T ] of the 2nd order simplex element s 146 146 A.3 [Q^] and [Ts] of the 3rd order simplex element 147 A.4 [Q^] and [Ts] of the 4th order simplex element 147 A.5 [Q^] and [T ] of the 5th order simplex element 148 A.6 149 s and [T ] of the 6th order simplex element s ix List of Figures 1.1 Three-phase cables and tunnel-installed single-core (SC) coaxial cables . 1 2.1 A (iJT+l)-conductor system 10 2.2 Circuit representation of a multiconductor system 10 2.3 An equivalent network for the line 21 2.4 A SC coaxial cable with K cylindrical conductors 24 3.1 Definition of simplex coordinates in a simplex with area S 30 3.2 The arrangement and location indices of the nodes in the fourth order simplex element (JV =4) 32 3.3 A commonly used node numbering scheme 33 3.4 [Q^] index string of the fourth order simplex element 34 3.5 Quadratic quadrilateral isoparametric element 37 3.6 Quadratic triangular isoparametric element 40 3.7 Matrix structure of the final equations for [Z] calculations 44 3.8 Programflowchart for [Z] calculations with FEM 46 3.9 Geometry of a SC coaxial cable and its FEM solution region p 47 3.10 Radial divisions for SC coaxial cables 48 3.11 Current density distribution in the SC coaxial cable at 6kHz 52 3.12 Current density distribution in the SC coaxial cable at 60kHz 53 3.13 Maximum errors in [R] and [L] at different $ 54 3.14 FEM meshes at 60kHz for different error limits 55 4.1 63 A shallowly buried SC coaxial cable x 4.2 A deeply buried SC coaxial cable and its FEM solution region 67 4.3 Earth return current and maximum errors in \Z\ for different r& 69 4.4 J distributions in the earth at different frequencies 72 4.5 Replacing a deeply buried SC cable with a current 4.6 Perturbation coefficient Cp as a function of |r /p | 78 4.7 FEM mesh at 6kHz for r =2S 83 4.8 Earth return current and maximum differences in [Z] for different r& . . . 4.9 E field at 6kHz given by the analytical solution e b filament e e 77 85 88 4.10 J distributions in the earth at 6kHz from the three approaches 88 4.11 The layout of a tunnel installed cable and the FEM mesh at 600kHz . . . 93 4.12 J distributions in the earth at 600kHz from the FEM 95 5.1 Direct capacitances for multiconductor systems 100 5.2 Direct capacitances under DC condition 102 5.3 Errors in [C] calculation of a SC coaxial cable for different span angles . 112 6.1 230kV three-phase cable systems 115 6.2 Meshes around the cables at 60Hz for the two systems 116 6.3 J distributions in the earth at 600kHz with the earth current being 1+j'O A120 6.4 A 230kV PT cable system 121 6.5 Meshes for the PT cables at 6kHz 122 6.6 \A\ distributions caused by the loop current at 60Hz and 6kHz 124 6.7 | J\ on the inner surface of the pipe caused by the loop current 124 6.8 A sector-shaped cable 129 6.9 \A\ distribution in the sector-shaped cable caused by the loop current . . 130 6.10 Potential distribution in the sector-shaped cable 131 6.11 A two-layer stranded conductor 133 xi 6.12 | J\ distribution in the stranded conductor B.l A current filament buried in the earth xii Acknowledgements I wish to express my gratitude to my parents, who inspired and supported me to undertake this venture. I wish to express my sincere appreciation to my supervisor, Dr. Hermann Dommel, for his excellent guidance, encouragement, readiness to help, and patience. I also wish to thank Dr. Martin Wedepohl for sharing his vast experience and offering many constructive suggestions during my thesis work. I would like to thank the Transmission Engineering Division of B. C. Hydro and Power Authority for providing some of the cable data used in this thesis. The financial support of the System Planning Division of B. C. Hydro and Power Authority and of the University of British Columbia is gratefully acknowledged. I would like to thank Dr. Jose Marti and my colleagues Luis Naredo and Liu Guang for many inspiring and helpful discussions, and to Dr. Luis Marti for letting me use his earth return impedance routine. Appreciation is also expressed to the staff of the Department of Electrical Engineering, UBC, for their assistance, especially to Mr. Robert Ross and Mr. Dave Gagne for their patience and help in answering so many of my questions. Last but not the least, I thank my wife, Jie, who through her incredible thoughtfulness, encouragement, and patience has helped me with every step during my doctoral program. xiii Dedicated m Jiy p a r s x m xiv Chapter 1 Introduction Underground power cables are widely used in electric power systems, and are therefore one of the important components in power system transient analysis. Unlike overhead power lines, cables are geometrically compact. The distances among conductors are comparable with the dimensions of the conductors. The conductors and cables often have non-coaxial geometries (Fig. 1.1), which make closed-form solutions difficult or impossible. Proximity effects are no longer negligible in cables, compared to overhead lines. Because of proximity and skin effects, cables possess much stronger frequencydependent characteristics than overhead lines, and mathematical models are needed for cables in the transient analysis which include these characteristics. pipe-type sector-shaped insulation sheath air earth core SC coaxial cables (a) Different types of three-phase cables (b) Tunnel-installed SC coaxial cables Figure 1.1: Three-phase cables and tunnel-installed single-core (SC) coaxial cables 1 Chapter 1. Introduction 2 Several accurate mathematical models for the transient analysis of power cables have been developed. For frequency domain based methods a general model for overhead lines and underground cables was developed by Wedepohl et aim 1969 [6]. As the transients are calculated in the frequency domain, the frequency-dependence can easily be included. For time domain transient analysis a general model was developed by J. Marti in 1982 [29]. In this model, curve-fitting techniques are used to find approximate rational functions to replace the frequency-dependent characteristics of mode parameters in the frequency domain. The mode parameters are related to the phase parameters through mode-phase transformation matrices. A similar model was developed by Humpage et aim 1980 [25]. In that approach, the z-transformation technique is used, and the curve-fitting is done in the z domain. These two models, however, do not consider the frequency-dependent characteristics of the mode-phase transformation matrix. J. Marti's model was improved by L. Marti in 1988 to take the frequency-dependence of the mode-phase transformation matrix into account [39]. The frequency-dependent phase parameters of the cable, i.e. series impedance matrix and shunt admittance matrix per unit length, are the input data for all the models mentioned above. The calculation of these parameters, however, is not an easy task, especially the calculation of the series impedance matrix. There are mainly two ways to calculate these parameters: analytically and numerically. In the early days before the arrival of computers, many theoretical studies were done on the calculation of series impedances of overhead lines and of underground power cables. In 1926 Carson derived formulas for the earth return impedance of overhead lines [1], and Pollaczek derived similar formulas for underground cables in the same year [2]. In 1934 Schelkunoff gave the complete solution for the electromagneticfieldsinside a single-core (SC) coaxial cable and derived corresponding impedance formulas[3]. The earth return impedance formulas and the impedance formulas for SC coaxial cables are still widely Chapter 1. Introduction 3 used in modern cable parameter calculation programs [10][23]. With the recent increase in computing power, more sophisticated formulas were developed as well. Tegopoulos et at derived the field solution for pipe-type (PT) cables with circular inner conductors in 1971 [8]. The formulas for the impedances of these PT cables were derived by Brown et al in 1975 [13]. The closed-form formulas, however, are limited to cables with simple geometries. Even for PT cables with circular inner conductors, several approximations have to be made in order to derive the impedance formulas. If irregular geometric shapes are involved, or if all the skin and proximity effects should be considered in the parameter calculation, then the closed-form formulas are no longer applicable, as there is no analytical solution for those problems. Numerical methods have to be used instead. In order to calculate the parameters of underground power cables with numerical methods, the corresponding electromagnetic field problem is solved first, and the parameters are then derived from the numericalfieldsolution. Based on the assumptions discussed in Chapter 2, the related electromagnetic fields are split into static electric fields, from which the shunt admittance matrix is calculated, and quasi-static magnetic fields, from which the series impedance matrix is calculated. Because of skin and proximity effects, it is more difficult to solve quasi-static magneticfieldsthan to solve static electric fields. Several numerical methods can be applied to solve the quasi-static magnetic fields needed for the series impedance calculation: subdivision method, finite element method (FEM), boundary element method (BEM), and hybrid method. The subdivision method divides conductors of a multiconductor system into small subconductors such that the current density within the subconductors can be assumed as uniform. A set of linear equations is established by using self and mutual impedance formulas among the subconductors. The series impedance matrix of the multiconductor Chapter 1. Introduction 4 system can then be found from the coefficient matrix of the linear equations by matrix reduction. The method wasfirstapplied to the series impedance calculation of multiconductor systems by Comellini et aim 1973 [9]. In their paper, round subconductors with uniform current density were used. Rectangular subconductors and subconductors with other shapes were used by Deeley et al [16] and Lucas et al [17] in 1978, as well as by Weeks et al for microstrip lines in 1979 [22]. This method was applied to the impedance calculation of PT cables by Arizon et aim 1988 [38]. With the subdivision method, the series impedance matrix of a multiconductor system is calculated without solving the field equations. Only the conductor regions are needed. The coefficient matrix, however, is a full matrix, and a homogeneous permeability in the space must be assumed. Because uniform current density in each subconductor is generally assumed, a veryfinesubdivision has to be used to achieve accurate results if strong skin and proximity effects are present. With thefiniteelement method, the domain of a closed boundary problem is divided into small elements such that the original unknownfielddistribution in each element can be approximated by certain functions expressed in terms of unknownfieldvariables at element vertices. Thefieldvariable values at all the vertices of the FEM mesh can then be determined with variational technique or with Galerkin technique. In the early 1970's, the method began to be applied to eddy-current related problems in power apparatus. It was used tofindthefielddistribution in magnetic structures such as motors by Chari in 1974 [12], in a case where super-conductivity was assumed for the current carrying conductors. The method was applied to study the skin effect in a single current carrying conductor by Chari et aim 1977 [14]. It was extended to study the skin effect in multiple current carrying conductor systems by Konrad in 1981 and in 1982 [26][28]. In his papers, source current density Js was related to measurable conductor currents; therefore, conductor currents replaced Js to become the forcing function. It was also suggested in his papers Chapter 1. Introduction 5 that the impedance matrix of a multiconductor system could be calculated from the J$ vector directly. In 1982 Weiss et al combined the original FEM equations, which had an unknown Js vector, with the equations which relate the Js vector to the conductor currents [31][32j. As a result, the field solution can be found in a single step. Thefinallinear equations in FEM generally have a symmetric, diagonally dominant, banded, complex coefficient matrix. It is easy for FEM to handle problems with regions having different permeabilities. The method is also flexible with respect to the shape of the elements and to the order of the approximating functions. Open boundary problems, however, cannot be handled easily by FEM, though the ballooning technique can be used for problems with Laplacian exterior regions. Also, the whole problem domain within the closed boundary has to be discretized in FEM. Instead of discretizing the whole problem domain, only the boundary of the problem domain is divided into small elements with the boundary element method, over which approximating functions with unknown coefficients are assumed. These unknown coefficients are found by applying Green's theorem and the basic solution for impulse sources (Green's function). There are several applications of BEM to eddy-current related problems [24] [33]. It has also been used to solve the field distribution in power cables [36]. However, few applications relate to the parameter calculation of multiconductor systems. BEM handles open boundary problems with Laplacian exterior regions easily. It has fewer variables because only the boundary is discretized. Thefinalcoefficient matrix, however, is generally an unsymmetric full matrix, and there is no general procedure to derive the Green's functions for arbitrary problems. It may be difficult to apply BEM to the earth impedance calculation where the cable systems are surrounded by the poorly conductive earth occupying half space. The hybrid method has been suggested to combine the advantages of FEM and BEM. Several applications to eddy-current problems have been reported [30][35]. Chapter 1. Introduction 6 In this thesis project, FEM is applied to the calculation of the series impedance matrix and shunt admittance matrix of underground power cables. Although FEM has been used to solve eddy-current related problems, most applications have been concerned with the field distribution in the whole problem domain at low frequencies. Only a few studies have been made on the parameter calculation of multiconductor systems[41] [42]. The following topics were studied: selection of element types; selection of function orders for high order approximating functions; optimal mesh generation; model development for infinite earth region in FEM; algorithm for solving the final banded complex matrix; and error analysis for different types of elements and for approximating functions with different orders. In Chapter 2 the principal equations for the calculation of the series impedance matrix ([Z]) of underground power cables are derived. The FEM based on the Galerkin technique is used to solve the principal equations. Two formulations for the [Z] calculation of a multiconductor system from the field solution, the Js method and the loss-energy method, are derived. These formulations are new in the literature, although the concepts already exist. In Chapter 3 two types of elements in [Z] calculations are discussed. They are high order simplex elements and quadratic isoparametric elements. The accuracy in [Z] calculations, the computational efficiency, and the mesh generation schemes of these elements are studied numerically by calculating the [Z] of SC coaxial cables. The Js method and the loss-energy method are also compared numerically. The integration matrices of simplex elements are given in the exact fraction form for the first order up to the sixth order. In Chapter 4 a new technique is developed to include the infinitely large and conductive earth in the FEM in [Z] calculations. This technique reduces the solution region drastically. Chapter 1. Introduction 7 In Chapter 5 the principal equations for the shunt admittance matrix ([Y]) of underground power cables, i.e. Laplace equations, are solved with the FEM, and the formulations to calculate [Y] from the field solution are discussed. In Chapter 6 the numerical results of several case studies are presented. These studies include the applications of Pollaczek's earth return impedance formula to multiphase underground cables with irregular structures; comparisons of parameters of PT cables or of sector-shaped cables between the results from the FEM and those from the analytical formulas; and the internal resistance calculation of stranded conductors with the FEM. In Chapter 7 the conclusions of this thesis project are presented and possibilities for future research work are discussed. Chapter 2 Impedance Calculation with Finite Element Method 2.1 Introduction This chapter is mainly concerned with the formulation of the field equations and the corresponding FEM solution for the series impedance of underground power cables. The line parameters, i.e. series impedance matrix per unit length ([Z]) and shunt admittance matrix per unit length ([i^]), are defined in connection with the transmission fine equations. The principal equations describing the magnetic fields for the [Z] calculation are derived. These equations are solved by FEM based on the Galerkin technique. In order to retrieve [Z] from the field solution, two methods are suggested. One is to relate the source current density vector [Js] to [Z] directly, and the other is to split the power loss and the stored magnetic energy in the system into summations of elements which relate to the elements in [Z]. As the magnetic fields in SC coaxial cables can be solved analytically, the analytical solutions can be used to check the accuracy of the numerical approaches. In this chapter, a general form of [Z] for SC coaxial cables is given, which will be used frequently in later chapters. 2.2 Definition of Transmission Line Parameters An overhead transmission line or an underground power cable can be represented as a general multiconductor system. The following assumptions are made for such a system. 8 9 Chapter 2. Impedance Calculation with Finite Element Method 1. The system is composed of infinitely long metallic conductors and the earth. The axes of the conductors are parallel to each other and to the surface of the earth. 2. The system is isotropic, linear, and longitudinally homogeneous. All the conductors and dielectrics have constant permittivity c , permeability /z , and conductivity a. P P So does the earth. 3. There is no volume charge inside the conductors and the earth. The charges are only located on the surfaces of the conductors and the earth. 4. Displacement currents in the conductors and the earth are ignored. 5. The frequencies used in the study are far below the value where the corresponding wave length A becomes comparable with the lateral dimensions of the system. With the above assumptions, a unique relationship betweenfieldquantities, E and H, and circuit quantities, V and I, can be established. Therefore, the electromagnetic fields in the system can be represented by a distributed-parameter electric network. The transmission line equations (telegraph equations), instead of Maxwell's equations, can then be used to describe the system. Fig. 2.1 shows a longitudinal section with infinitesimal length dz of a (.K"+l)-conductor system. The (K+l)th conductor is used as a reference conductor. V\, V , ..., VR are 2 the conductor voltages with respect to the reference conductor K+l, and / j , 7 , . . . , IK 2 are the conductor currents. They are phasors. The z axis is in parallel with the axes of conductors. The circuit representation of the system in the frequency domain is shown in Fig. 2.2. In Fig. 2.2, all the variables are in vector form, and all the parameters are in matrix form and are quantities per unit length. [i2(u>)] and [Lc(^)] account for the power loss and magnetic energy storage in the conductors, respectively. [LD] accounts for the magnetic Chapter 2. Impedance Calculation with Finite Element Method 10 energy stored in the dielectrics. [C] and [(?] represent the electric energy storage and the power loss in the dielectrics, respectively. conductor L _ - - J g ^ , conductor k . {*+«'* Figure 2.1: A (iif+l)-conductor system [/] [L )dz [R((a)]dz {1^(0!)] dz D O [V] [V]+d[V] [C]dz [C\dz I Figure 2.2: Circuit representation of a multiconductor system The transmission line equations corresponding to Fig. 2.2 are d[V] = dz ([R(u,)}+ML(u)])[I] = [Z(u)][I\ (2.1) M dz ({G]+MC])[V) = [Y(u)}[V} (2.2) = in which [V] =[V V ,...,V ? U 2 K (2.3) 11 Chapter 2. Impedance Calculation with Finite Element Method [I] =[I I ,...,I ] U 2 (2.4) T K [Z(w)] = [R(v)] + ML(u)] = [J2(w)] + M[L (u>)] + [L ]) c [*»] =[G]+MC] D (2.5) (2-6) [Z(u;)] and [Y"(w)] are the series impedance matrix per unit length and shunt admittance matrix per unit length of the system, respectively. They are generally referred to as transmission line parameters. In most cases, [Z(u>)] is a nonlinear function of u>, while [Y"(o;)] can be simplified as ju>[C] by ignoring [G]. For simplicity, u will be omitted from [Z(u)] and [Y(u>)]. The transmission line equations are the fundamental equations on which all transient analysis models for underground power cables are based [10] [29] [39]. [Z] and [V] are the input data for those models. According to the assumptions made before, [Z] and [Y] can be calculated from a quasi-static magneticfieldand a static electricfieldof the system, respectively. 2.3 Principal Equations in Impedance Calculations The matrix [Z] of an underground power cable can be found from the magnetic field distribution in and around the cable. In order to do so, the principal equations describing the magneticfieldshave to be derived first. Based on the assumptions 1 and 2 given in the previous section, the electromagnetic field of a cable system is two-dimensional and linear. In order to isolate the magnetic field of the system, displacement currents in dielectrics are also ignored in addition to assumption 4. As a result, the original electromagneticfieldbecomes a quasi-static magnetic field which is excited solely by the conductor currents. Ignoring displacement currents in the series impedance calculation is common practice[l][2]. To justify this assumption, Wedepohl and Efthymiadis rigorously analyzed the overhead line case over the full frequency 12 Chapter 2. Impedance Calculation with Finite Element Method spectrum[18][19]. They found that the assumption was valid until the height of the conductor above ground became comparable with 1/4 of the wavelength. The same result will hold for a cable system, but since the dielectric spacings are only a few centimeters, differences will only occur above several hundred MHz, which is orders of magnitude above the frequencies of interest in power system transients. Recent comparisons betweenfieldtests and simulations also confirm the validity of the assumption[39][43]. The displacement currents are therefore ignored in the series impedance calculation. With the discussed assumptions, the following equations are derived from Maxwell's equations V x E = -JUJB (2.7) —V x B = J (2.8) V E =0 (2.9) Introducing the magnetic vector potential A as B = VxA (2.10) and inserting (2.10) into (2.8) and employing unity V x V x A = V(V • A) — V A , the 2 following equation is obtained - V ( V • A) - - V A = J (2.11) V A =0 (2.12) - iv A =J (2.13) 2 ft ft Assuming (2.11) gives 2 As the current density has only a longitudinal component, J , E , and A can be written respectively as J u , Eu , and AVL . U is the unit vector along the z axis. By inserting z z Z z Chapter 2. Impedance Calculation with Finite Element Method 13 (2.10) into (2.7), the following equation can be derived (E + jwA)u = -V<f> (2.14) z where <f> is a scalar function. Reference [28] gives the following properties of <f> in a conductor: the equal value surfaces of <j> are perpendicular to the z axis and V<f> is a constant within each conductor. These properties make it possible to define a unique voltage between one conductor and the reference conductor at a location along the system. Constant — V<£ in each conductor is defined as the source electric field Es. The physical meaning of Es (or —V<f>) is the voltage drop along a unit length of the system. The current density corresponding to Es is called the source current density Js, which is a constant over the cross section of a conductor. Three quantities are related by J u = (rE u = —<rVd> s z s z (2.15) (2.14) can now be written as J = -jvcrA + J (2.16) s Combining (2.13) with (2.16), the following linear two-dimensional diffusion equation can be derived - V A - jwvA + J = 0 (2.17) 2 s The integration of (2.16) over the cross section of a conductor gives 1=1 Jds = -ju f a Ads + S Js JSc c (2.18) JSc in which Sc is the cross-section area of the conductor. For a multiconductor system, there is one such equation for each conductor. (2.17) and (2.18) are the principal equations describing the quasi-static magnetic fields needed for the [Z] calculation. These equations are solved with FEM in the next section. 14 Chapter 2. Impedance Calculation with Finite Element Method 2.4 F E M Based on the Galerkin Technique for the Principal Equations For a multiconductor system with K conductors, the principal equations in the [Z] calculation can be summarized as + JS -V A-ju(rA 2 A* -Jul = 0 <7Ads + S Js Ck Js inS = h k (2.19) R {k = l,2,...,K) (2.20) Ck with boundary conditions ^lr 0 dA dn ri = 9o(x,y) (2.21) = 0 (2.22) SR is the solution region surrounded by boundary r = ro + IV Sc and Js are the crossk k section area and the source current density of the kth. conductor, respectively. go(x,y) is a known function. To and Ti are the Dirichlet boundary and the homogeneous Neumann boundary, respectively. 2.4.1 The Galerkin technique for solving differential equations The Galerkin technique approximates thefielddistribution A satisfying the differential equation (2.19) by a finite set of base functions ip (n = 1,2,..., N) as n N A = ip + Y, °n<p 0 where a , a ,..., ajv x 2 & i e (2.23) n unknown coefficients. ip is such that Q Mv =% (2.24) 0 V'lj ¥>2> • • -i VJV form a subset of the complete base functions of A. ip satisfies the n boundary condition Vn|r =0 o (n = l,2,...,JV) (2.25) 15 Chapter 2. Impedance Calculation with Finite Element Method Putting the approximate solution (2.23) into (2.19), a residual is introduced as R(A) = -V A - juaA + J /* in S 2 s (2.26) R The Galerkin technique forces this residual to satisfy the following integral equation / R(A)ip ds = 0 (n = 1,2,...,iV) n Js (2.27) R The above equation means that the projection of the residual on each base function is zero, or the residual is orthogonal to all the base functions. As there are N integral equations in (2.27), N unknown coefficients in (2.23) can be determined completely by (2.27). The forms of the base functions (also called trial or shape functions) will be discussed in the next chapter. In general, they should be simple functions. Polynomials are among popular base functions because they can be easily differentiated and integrated. With (2.23) and (2.26), (2.27) becomes / i (V • V(V>o + £ s J l*\ R £i a <p ) \ <p ds f J Js n n m jW(^ + £ 0 a ip )<p ds + n n m R + / Js Js<Pmds = Q (m = l,2,...,AT) (2.28) R Applying Green's formula V • (vVu) = Vu • Vu + vV u, the first term in the above 2 equation becomes r - / 1 r 1 -Vy? • V(V>o + £ a <p )ds + - V • (y> V(^ + £ nVn))ds N N a m S fi n J R n = = - / R^ JS i n m 0 JSRH N = 1 I^ -V^-V^o+Ean^+Z n=i °+ l Jr r I 1 ^ m ( 9 + a £ n=l (m = 1,2,..., AT) a J ^ r 9 n (2.29) It is shown in [44] that boundary conditions in (2.22) and (2.25) can be automatically satisfied if the second loop integral in the above equation is set to zero. Therefore, once <p n Chapter 2. Impedance Calculation with Finite Element Method 16 is established, the N unknown coefficients are found by solving the following equations (m = l,2,...,iV) (2.30) The differentiation order in the above equations is one order lower than that in (2.28). It is very difficult to use simple form base functions to approximate the original complicated field function over the whole solution region directly. If the solution region is divided into small subregions such that the original field function changes smoothly in each subregion, it will be possible to approximate the function by simple base functions, such as polynomials, within each subregion. This is the main idea behind finite element methods. 2.4.2 The F E M based on the Galerkin technique With the FEM, the solution region is divided into elements (subregions) and the base functions are systematically established in each element, as required by the Galerkin technique or other techniques for solving the differential equations. The algebraic equations (2.30) can be assembled element by element, and the unknown variables can then be found. The process for dividing the solution region is called meshing process, and the resulting region made up of elements is called a finite element mesh. All the element vertices and additional locations either on the element sides or inside the elements are defined as nodes. With the Galerkin technique, the values of the field variables at the nodes become the unknowns in (2.23). In the [Z] calculation, the field variable is A, and its nodal values are A (n = 1,2,..., N), with N being the number of nodes within the mesh. The n coordinates of A are represented by (x , y ). The shape function y> (n = 1,2,..., N) n n n in the [Z] calculation satisfies the following conditions. n 17 Chapter 2. Impedance Calculation with Finite Element Method 1. ip is continuous inside SR. Therefore, it is continuous across the interelement n boundary; 2. { 1 m=n (n,m = 1,2,. ..,N) 0 m^n When the solution region is divided, the boundary is also discretized. The boundary will be made of boundary element sides, and the nodes on the boundary element sides become boundary nodes. The original continuous Dirichlet boundary condition is now represented by those boundary nodes on r with prescribed values. A similar expression 0 can be written for ipo as N ^0 = YL ABi lPBi t=l B in which AB% (2.31) represents the known node value on the boundary To, <fiBi represents the corresponding shape function which has the same characteristics as <p discussed above, n and NB is the number of Dirichlet boundary nodes on T . 0 The nodes can be renumbered such that the first N nodes are the nodes with unknown node values and that the last NB nodes are Dirichlet boundary nodes with prescribed node values. The total node number is NT = N + NB. The approximate solution for A can be written as (2.32) A = £A < Y N PN n=l and (2.30) and (2.20) become r 1 r r - / -V<^ • V J* nVnds - I j<jJ<T(p J2 A ip ds + / J <p ds Js fi ^ Js ^ Js N t N t A m m R = - V) A n R N I n s m R (-V<p • V<p + ju}crip (p^\ ds+ I JsVmds = 0 m n m (m = l,2,...,iV) (2.33) Chapter 2. Impedance Calculation with Finite Element Method - J ^ M I n=l <7<Pnds) + S Js Ck k =h 18 (k = 1,2,..., A") (2.34) "C /S fc Equations (2.33) and (2.34) are combined together to form the following final algebraic equations in matrix form [U]+ju[T] [0] (2.35) \ -MGo] [ [F] , [F Y ] [So] T B in which U = ISR ^V<^ • V<p ds m = 1,2,....W T — fs VWmVnds n= M K = I FB* = Is fids [Gc] = diagfa,a ,... [Sc] = diaglScnSct^.-yScK] [A] = [Js] = mn F mn m n R SCk 1,2,...,N l = N + l,N + fmds T 2,...,N T Jb = l,2,...,JT Ck K [A A ,...,A ) U (2.36) ,cr ] 2 2 T NT Psi > Js ) • • • j ^5 ] r K 2 In matrices [F] and [F^] only the elements corresponding to the nodes in the conductor regions have non-zero values, and the others are zero. [Sc] is the conductor cross-section area matrix. [0] is a zero vector. [Gc] is the conductivity matrix. It is assumed that each conductor has uniform conductivity in its cross section; otherwise, [Gc] could not be extracted explicitly as in submatrix [Gc] [ [F] , [FB] ]• In practice, (2.35) is assembled T T element by element. That means all the integrals in (2.36) become the summations of integrals over the elements as [ JSR [ JS ( 2 37 ) Ei M is the number of elements in the finite element mesh. SE; (i = 1,2,...,M) is the region in the ith element. SR = SEI. U S&Z U • • • U SE , and Sc* £ SR. Inside the ith M Chapter 2. Impedance Calculation with Finite Element Method 19 element (2.32) becomes N EI A=E A*<p« in S (2.38) Ei n=l in which N Ei is the number of nodes in the ith. element, A ip is the shape function defined locally in the element. U Ei Ei U mn and T m n become • V(p ds = Ei is the node value of A, and Ei rpEi _ mn (2.39) ds JS Bi J? i — E for conductor elements Js Ei for conductor elements *%> = f <fif ds f where m = 1,2,...,N Ei excluding boundary nodes, n = 1,2,..., N , I = 1,2,..., A^;. Ei excluding unknown nodes, and i = 1,2,..., M. Boundary nodes are those on the boundary T with prescribed node values whose global numbers are bigger than N, and unknown 0 nodes are the nodes with unknown node values whose global numbers are equal to or smaller than N. The unknown nodes are made up of nodes inside the solution region and of those on the boundary IV fi and o~ are the permeability and conductivity in the Ei Ei ith element, respectively. The above formulas use local node numbers, while each node has a global node number. A node on an element side may be shared by several adjacent elements. When the assembly (2.37) is done, those elements will have contributions to the entry of the same global node in the final matrix (2.35). In (2.35), the unknown variables are the unknown node values and the source current density Js in all the conductors, while the forcing factors are the Dirichlet boundary node values and the conductor currents. Once the boundary conditions and the conductor currents are given, the magneticfieldrepresented by the discrete node values of A can be found by solving (2.35). From thefieldsolutions, the series impedance can be calculated. This is the topic of the next section. Chapter 2. Impedance Calculation with Finite Element Method 2.5 20 [Z] Calculations from the Field Solutions Two methods have been developed to calculate [Z] from the field solutions: the Js method and the loss-energy method. With the Js method, [Z] is calculated directly from the vector [Js]- With the loss-energy method, the power losses and the stored magnetic energy are calculated from thefieldsolution first. [R] is then calculated from the power losses, and [L] from the stored magnetic energy. The basic concepts on which these two methods are based have existed in the literature. However, they have not been applied to the impedance calculation of multiconductor transmission line systems in the way described here. The loss-energy method used in [41] for the impedance calculation differs from the method of Section 2.5.2. 2.5.1 The J method s As mentioned in Section 2.3, the physical meaning of Es in (2.15) from thefieldanalysis is the voltage drop per unit length along the system[28]. This corresponds to in (2.1) from the circuit analysis formulation. The relationship is fe] = - f ! (2.40) where [E ] = [E i,E , • • •, E ] - Therefore, T s S S2 SK [Js] = [G ][E ] = -[Gc]^C 8 If [/] is assumed as Ii = 0 and Ij ^ 0 = [Gc][Z][I] (2.41) = 1,2,..., K; i ^ j), the jth column of [Z] can be derived from the above equation as Zii = T i (« = 1,2,...,J0 (2.42) lj(Ti In order to find the whole [Z] matrix of a multiconductor system, the corresponding equation (2.35) has to be solved K times. Each time only one conductor carries current, Chapter 2. Impedance Calculation with Finite Element Method 21 while the other conductors are open-circuited. This method is simple and straightforward. There is no need to know all the field distributions before calculating [Z]. 2.5.2 The loss-energy method When the field distributions represented by A and J are known, the power loss and the stored magnetic energy in the system can be calculated. They are related to the series resistances and inductances in the system. However, as the system is made up of multiple conductors, the power loss and the magnetic energy of the whole system cannot be used directly in the [Z] calculation. In the following, the formulas for calculating [Z] from the power loss and the stored magnetic energy are derived by comparing the corresponding formulas in the circuit analysis and the field analysis formulations. In the derivation, the loss and the energy are broken down into summations of elements. These elements are then related to the elements in [Z]. Equation (2.1) can be represented by the equivalent network shown in Fig. 2.3. For such a resistive-inductive passive network, the complex ~ ST Figure 2.3: An equivalent network for the line power going into the network is S = P-jQ = Chapter 2. Impedance Calculation with Finite Element Method 22 where S, P, and Q are time-average complex power, time-average power loss, and timeaverage reactive power in the equivalent network. Assuming [I] = [IR] + j [2j] and using the fact that P = P and Q = Q , the following equations can be derived T T P = = T T PV [IR] [R}[IR] + [II] [R][II} T T i t=l j=l K Q (2-44) I K = o,([//[i][/ ] + [/;f[I][//]) = E E ^ (2-45) fl t=i i=i in which Pij = RiiWRi + hh) = UL^IRJR, Q I I + (2-46) //.//.) (2.47) The time-average magnetic energy WM stored in the system is related to Q by WM = ±Q = hlR} [L][I ] + T R [//] [I][//]) = E E w T Mii (2.48) in which wiiii^LiiWRi+hhi) Consequently, if (2-49) and WM^ can be calculated from the field solution, Rij and Lij will be given by (2.46) and (2.49), respectively. The formulas to calculate the time-average power loss and the time-average magnetic energy stored in the fields are ' K J J* ~ir p = Y.L fc=l W M ( - °) ds JS c k 2 5 <* = I I BH*ds = lT Re( 2Js / R 2£i [ \Jsc k AJ*ds\ J (2.51) where all the phasors are in RMS. These formulas, however, only give the total power loss or the magnetic energy in the system. When the system has an arbitrary combination of Chapter 2. Impedance Calculation with Finite Element Method 23 conductor currents, J in a conductor is caused by all the conductor currents. Because the system is linear, J in the conductor can be written as a summation of different current densities caused by different conductor currents. Assuming J( ) is the current density in K the kth. conductor, it can be written as K (A: = 1,2,..., it) J(K)=T, (KI) J t=i where (2.52) is the current density in the kth. conductor caused by the current /,• in the ith. conductor. (2.50) now becomes k=l c K . Js k=\ c a Js k = E / k=l Ck K K JS -i K ; s a k K E E K ,=i K , W=E E E / 3<=1 j=\ a i=i K 7.,.. ™ S i=l 7=1 k=l C„ JS * a = EE ^J (2-53) »=ii=i = — E J* / P I J ™ k=l k k=l c ^ (2.54) Jbc Js a k Similarly, = £ £ o £ .'=1 j=l ^ where = ^E = J2iLi ^.(Ai) Z R e R U = E E ^ V- ^ / 76 A=l (^ e «=1 j = l V-) (V ) J c (- ) 5 2 Kk=l k La = ^ds)l{I I +I I ) Jjt JSc a 56 is the A in the Mh conductor caused by current /,-. Relating (2.54) and (2.56) with (2.46) and (2.49), respectively, gives **i = {EL (2.55) J JU Ri Ii I} and L{ - as } {i,j = l,2,...,K) ^«)^y)*)}/( ^ *i+V/i) J J (2.57) (i,i = l,2,...,A')(2.58) 24 Chapter 2. Impedance Calculation with Finite Element Method In order to obtain J(^,) and (k,i = 1,2,..., K), (2.35) of the system has also to be solved K times with only one conductor carrying current at a time. Compared with the Js method, the loss-energy method needs the full information about the field distributions at least inside the conductors in order to retrieve the [Z] matrix. It also requires additional calculations to evaluate the integrals in (2.57) and (2.58). The two methods are applied to simple coaxial conductor systems in Chapter 3, and the comparisons are made for the elements of [Z] with different kinds of elements. 2.6 Impedances of Single-Core Coaxial Cables The quasi-static magnetic fields inside single-core (SC) coaxial cables can be solved analytically due to axisymmetrical geometry of the cables. The field solutions were derived by Schelkunoff in 1934[3]. The formulas for the impedance calculation of SC coaxial cables can then be derived from the analytical field solutions. In later chapters, these impedance formulas are used frequently to test various numerical approaches. Fig. 2.4 shows a SC coaxial cable with K cylindrical conductors. r^k and rsk are internal and external radii of thefcthcylindrical conductor, respectively. With the Chapter 2. Impedance Calculation with Finite Element Method 25 same assumptions as made in Section 2.2, the transverse magnetic field Hug and the longitudinal electricfieldEu satisfy the following equations in cylindrical coordinates z E =L4±m ( 2 . 6 0 ) ro~ dr The solution for E inside the Arth conductor is E{r) = — {C h(r/Pk) - C K (r/p )) o~Pk Ik Kk 0 r k Ak <r<r (2.61) Bk where C Ik = 7T7^{—Ki(rWP*) + —KxCrxJp*)) ZTrL>£) \r k CKH = « rB Ak ^ f^Ii(rB*/w) + ^IiCrxik/p*)) ZnCDk \rAk CDk = Pk I Ak and I a r e Bk ) k r (2-62) ) Bk h(rBk/Pk)Ki(r k/Pk)-IiirAk/P^Kxirek/pk) A = \ -r^ V JUCTkHk internal and external return currents of thefcthconductor, respectively. p is a complex penetration depth in the kth conductor. I , and K (n = 0,1) are the k n n modified Bessel functions of the nth order and respectively of thefirstand second kinds. E on the internal and external surfaces of thefcthconductor can be derived as E(r ) = Z I Ak E(r ) Ak Ak Ak +Z I Mk =Z I Bk Mk (2.63) Bk +Z I Bk (2.64) Bk where Ak = 2TrrAkl p C Dk Bk = Dk z z k k k k 1 + (^( + I 2irr a p C Bk ( °( '^ /^) ( fe/P ) 7 fe rBfe Kl rfl fc Kl r fe (2.65) o( Bk/Pk)h(r /p )) (2.66) r B /P ) ( ^/P ) fc o( Ak/Pk)li(r k/Pk)) K K r Ak k , 26 Chapter 2. Impedance Calculation with Finite Element Method ZAICI Zsk* a n < i %Mk defined by Schelkunoff as internal surface impedance, external a r e surface impedance, and transfer impedance of thefcthcylindrical conductor, respectively. Suppose that the system in Fig. 2.4 has a superconductive surface at TA . K+1 This surface will be used as the reference for the conductor voltages and the return path for the conductor currents. The series impedance matrix of the system can then be derived as [10] [23] t Zd yod 7. 1 L yod ^2 L Z3od [Z] \ L yd 1 yod "3 z$ yod 3 z$ L yod ^3 yod K L \ (2.68) z$ yod K yod K L Zjc J where K zt = K Yl EQ Z K - 2 E k-i K zf J K = = %Bk + %Dk + %Ak+l Z ZBHI BK and Z\ik (i = l,2,...,K) (2.69) K 2TT ZAk+ii Mk k=i+l ZzQ-k ~ ZMX - 2 ^] Z\fk EQk ZEQ Z a r e +Z \ (i = 2,3,..., K) (k = l , 2 , . . . , t f - l ) (2.70) (2.71) (2.72) 2 r Bk , (i = l,2,...,K) (2.73) given by (2.65)-(2.67). Zm is the impedance related to the inductance in the insulation between the kth. conductor and the (k -f l)th conductor. 2.7 Summary In this chapter the basic assumptions in the impedance calculation of underground power cables with FEM are discussed, and the principal equations are derived from Maxwell's Chapter 2. Impedance Calculation with Finite Element Method 27 equations. The principal equations are then solved with FEM using the Galerkin technique. The Js method and the loss-energy method are derived for the impedance calculation from the field solution. The general form of [Z] for SC coaxial cables is also given in this chapter. Chapter 3 Element Types and Shape Functions 3.1 Introduction The finite element formulation derived in the last chapter is incomplete because the shape function <f> ' (n = 1,2,..., A T E J inside an element is still unknown. The shape E function <f> is related to element shapes and interpolation functions in the elements. In Ei this chapter two kinds of elements are discussed in connection with [Z] calculations: the simplex element and the isoparametric element. The corresponding shape functions of these elements are used to complete the finite element formulation. In most two-dimensional eddy-current related problems, high-order simplex elements are used to calculate the magnetic fields. These problems are generally dealing with the fields in non-current-carrying conductor regions or non-conductor regions at low frequencies. In [Z] calculations, however, the frequency will vary from 0 to 1 MHz, and accurate solutions for the fields inside the current-carrying conductor regions are very important to obtain [Z] with good accuracy. At high frequencies, the conductor currents are concentrated in narrow regions near the conductor surfaces, which generally follow the contours of conductor surfaces. Under this situation, the curve-sided isoparametric element may be more suitable than the straight-sided simplex element due to the fact that most underground power cables are made from round conductors. Very little work has been done so far on the shape functions for [Z] calculations. To 28 Chapter 3. Element Types and Shape Functions 29 find out which, element is more suitable, the high-order simplex element or the curvesided isoparametric element, several factors have to be considered: the accuracy of [Z]\ the computation efficiency, which depends on the number of nodes in the mesh, on the final matrix bandwidth and pivoting, and on the CPU requirement; and the mesh generation. In this chapter, SC coaxial cables are used to study the shape functions in [Z] calculations. Because the [Z] for a SC coaxial cable is known, the FEM results with different shape functions can be compared with the theoretical results. This is used as an accuracy criterion. The factors mentioned above are compared for different elements. Also, based on the numerical results, the J$ method and the loss-energy method for calculating [Z] from thefieldsolutions are compared. 3.2 High-Order Simplex Elements The approximatefieldsolution given by (2.32) in the preceding chapter can be viewed as a two-dimensional interpolation formula. Once thefieldvalues at discrete nodes are known, thefieldvalue at an arbitrary location within the solution region can be calculated with the interpolation formula. The shape functions in (2.32) can also be established through setting up such an interpolation formula. As briefly mentioned in Section 2.4.1, the polynomials are among the popular shape functions used infiniteelement analysis, because they can be differentiated and integrated easily. The corresponding shape functions of the simplex elements and the isoparametric elements discussed in this section and in the next section are based on polynomials. Therefore, the shape functions to be discussed can be treated as interpolation polynomials which satisfy the conditions given in Section 2.4.2: the shape function for the [Z] calculation must be continuous across the interelement boundary; it must reach unity at 30 Chapter 3. Element Types and Shape Functions its associated node and it must be zero at other nodes. A simplex in N-dimensional space is the minimal possible nontrivial geometric figure defined by N +1 vertices. Therefore, a simplex in two dimensions is a triangle. A triangle (simplex) with area S is shown in Fig. 3.1(a). An arbitrary point P inside the triangle can be used to split the triangle into three subtriangles with areas Si, S2, and 53 as shown in thefigure.Simplex coordinates are defined as C< =| * = 1,2,3 (3.1) The location of point P can be uniquely defined by the simplex coordinates. • x (a) global coordinates (b) simplex (local) coordinates Figure 3.1: Definition of simplex coordinates in a simplex with area S From (3.1), the simplex coordinates are obviously varying in the range between zero and one. They are independent of the location of the triangle in the original x-y coordinate system. Only two of the three simplex coordinates are independent because Ci+e +c3 = 2 5 l + ^ + 5 3 = i (3- ) 2 With (3.1) a triangle in x-y (global) coordinates in Fig. 3.1(a) will be transformed into another one in simplex (local) coordinates in Fig. 3.1(b). 31 Chapter 3. Element Types and Shape Functions The interpolation polynomials required by FEM for [Z] calculations can be easily established with simplex coordinates. The detailed derivations can be found in [7] and [34]. Some key equations are given below, and the integrals in (2.39) are completed. 3.2.1 Shape functions in simplex elements Applying the formula calculating the area of a triangle from the coordinates of its vertices to (3.1), the relationship between the simplex coordinates, or local coordinates, and x-y coordinates, or global coordinates, can be found as ( 1 1 = 25 ) 322/3 - xy 2 2/2 - 2 / 3 «3 - «2 332/1 - «lJ/3 2/3 - 2 / i x - x ^ Zll/2 - Z2J/1 2/i - 2 / 2 x 3 x 2 - \ X 3 3i / \y ) In a simplex the shape functions can be expressed in terms of simplex coordinates as [34] m +m +m (N ,Ci)P (NP,C2)Pm (NpiC3) m2 p in which P (N ,Ci), mi 3 and P P (N ,C ), m2 p P (N ,C) m p p p 3 = AllW-fc) m ! Po(AT ,C) = p (N ,£ ) m3 2 1 2 = 3 N p (3.4) are auxiliary polynomials given by m = l,2,...,AT p (3.5) fc=0 1 (3.6) N is the order of the polynomials. Integer indices n»i, m , and m have values from p 2 3 0 to N and are related to the locations of nodes inside the triangle. They will be p called location indices. The nodes inside a simplex element are at the intersections of three groups of equally spaced parallel lines as shown in Fig. 3.2(a). Each group is in parallel with one of the triangle sides. The lines in a group are numbered from 0 to N p starting from the line aligned with the triangle side. Therefore, for each node there are three numbers associated with three intersecting lines where the node is located. These 32 Chapter 3. Element Types and Shape Functions numbers indicate the locations of the lines in different groups and are the location indices minima in (3.4). The location indices can be easily extracted from Fig. 3.2(a) into an explicit form as shown in Fig. 3.2(b). (a) element node arrangement (b) the location indices of the nodes Figure 3.2: The arrangement and location indices of the nodes in the fourth order simplex element (iV =4) p 3.2.2 Integral matrices of simplex elements In practice, the nodes in a simplex element can be numbered in any order; however, each node will have the fixed location indices. A commonly used node numbering scheme is shown in Fig. 3.3. For this scheme, the shape functions in (2.38) are related to the shape functions in (3.4) as y?f* = ct^, ip = a( _!) ,..., and ip^ = Ei n 10 aon' For the nth order 0 polynomials, NE = (n + l)(n + 2)/2. { With such a node numbering scheme and the shape functions given by (3.4), the integral UJH\ and T% in (2.39) can be derived as[34] n = 7-E (3-7) 33 Chapter 3. Element Types and Shape Functions N -l El Figure 3.3: A commonly used node numbering scheme (3.8) in which V m n Js* \dC dCk-x) \d( i k+l JS , E< (k = 1,2,3) dC -iJ 2S k (3.9) Ei (3.10) &E: A 0 k+ is the included angle of vertex k in element E{. As all the shape functions are in the form of (3.4), the integration functions in (3.9) and (3.10) are the functions of simplex coordinates only. They are the same for all the elements, independent of the element shapes and sizes. Therefore, they can be evaluated once and for all and tabulated. From (3.3), the following equation can be derived dx d3 = dxdy = 8x By By S<1 B< d(id£ 2 (3.11) = 2SEidCid£ 2 2 As the integral functions in (3.9) and (3.10) can be broken into summations of product terms of simplex coordinates, these integral equations are evaluated by ds r 1 r -^ 1 ...... + * + 2)! (3.12) 34 Chapter 3. Element Types and Shape Functions Under the node numbering scheme shown in Fig. 3.3, the matrices [ Q ^ ] and [Ts] of the first order to the fourth order have been given in the exact rational form in [7]. These matrices and those of the fifth and the sixth orders are listed in Appendix A as a reference. There are two typesetting errors in [7], one in [QW] of the third order and one in [Ts] of the fourth order. [QW] and [QW] are related to [QM] by = Qo}m)0(n) Qmli = Qo(m)0(n) = Qolo(m))0{0(n)) 0() is the index string for [0, W]. Fig. 3.4 shows how to find the index string for the fourth order simplex element. The vertices have to be numbered counterclockwise. The nodes ® Figure 3.4: [ Q ^ ] index string of the fourth order simplex element in the simplex are numbered in the scheme shown in Fig. 3.3 starting from vertex one. The numbers are those without parentheses in Fig. 3.4. The nodes are renumbered in the same way starting from vertex two. The numbers are those with parentheses. Then the first group of numbers are to be the indices in 0(), and the second group of numbers are to be the corresponding values of 0(). For the fourth order, 0(1) = 15, 0(2) = 10, 0(3) = 14, and so on. Tab.3.1 shows [ Q ^ ] index strings for different orders. Chapter 3. Element Types and Shape Functions 35 Table 3.1: [QW] index strings Index string 0() Order l2 3 6 10 15 21 1 28 3 3 4 5 6 1 3 6 10 15 2 21 5 2 5 9 14 20 4 27 8 1 3 6 10 7 15 12 2 5 9 14 11 20 17 4 8 13 19 16 26 23 1 3 6 2 5 9 4 8 13 7 12 18 1 3 2 5 4 8 7 12 11 17 10 1 14 2 19 4 25 7 6 11 9 16 13 22 18 24 F^ and F§* can be calculated from [Ts] as k k *S = S J2Ts Ei (3.13) mn n=l = (- ) 3 14 n=l where m = 1,2,... ,NE , excluding boundary nodes, and I = 1,2,... ,NE , excluding T { unknown nodes. Obviously, they can be calculated from the tabulated [Ts] directly. No integral calculation is needed. High order simplex elements satisfy the condition of interelement boundary continuity. With the nth order shape function, the approximation in (2.38) will also be an nth order polynomial which is uniquely defined by (n + 1) nodes on each element side. Therefore, if two adjacent elements have the same node value for each node on the shared element side, the function A will be continuous across the side. High order simplex elements are simple. The integrations are exact and independent of triangle shapes and sizes and can be done once and for all. The nodes for higher orders can be easily created from the first order element mesh by the computer. Required by the shape function, the element sides are straight, and the nodes are evenly located. This may not be suitable for [Z] calculations of underground power cables because most cables are made of round conductors. In the next section, the isoparametric element with Chapter 3. Element Types and Shape Functions 36 curved sides will be discussed. 3.3 Isoparametric Elements Cylindrical shells or round solids are common conductor shapes with underground power cables. Although elements with straight sides such as simplex elements can always be used to mesh a region with curved boundaries, many elements are required to fit the mesh into those boundaries. If elements with curved sides are used, instead, the number of elements in the mesh can be reduced. For simplex elements, it can be seen from (3.3) that the relationship between the global coordinates x-y and the simplex (local) coordinates C1-C2-C3 is linear. Therefore, a straight side element in the local coordinates will be mapped into a straight side element in the global coordinates. If a non-linear relationship is used, a straight side element in the local coordinates can be mapped into a curved side element in the global coordinates. In this section, two types of curved side elements are briefly discussed: the quadratic quadrilateral isoparametric element and the quadratic triangular isoparametric element. Detailed derivations can be found in the literature [27] [37]. 3.3.1 Quadratic quadrilateral isoparametric element In an isoparametric element, the global coordinates x and y are expressed as the interpolation functions of x coordinates and y coordinates, respectively, of the element nodes. These interpolation functions are in the same form as those used for the field variable in the element. In Fig. 3.5, a quadratic quadrilateral element is shown in both global and local coordinates. If the element nodes are numbered as shown in Fig. 3.5(a), interpolation formulas for A, x, and y in element Ei can be written as A = (3.15) 37 Chapter 3. Element Types and Shape Functions © f (-1.D ^ ^ D , . (0.1) E © '© © © (1.0) (-1.0) © © © (-1.-D (1.1) (0.-1) 1) © (1,-1) (b) local coordinates (a) global coordinates Figure 3.5: Quadratic quadrilateral isoparametric element x = \j3][* ) (3.16) y = \P}[y ) (3.17) Ei Ei where [A <] E [**] — F l >•2 > • far*] = [yfSyfv c ..,*?]* (3.18) [/3] is a row vector, and its elements are denned by the following expressions in the local coordinates [37]. /?! = J ( l - W ) ( l - I / ) ( ~ 1 - V - I / ) f3 = \{l + 2 & = 1(1 - v )(l - *) 2 v)(l-v)(-l+v-v) (3.19) P = i(l + v)(l + u)(-l + v + u) fh = i(l-v'){l P = \(l-v)(l /? = 1 ( 1 - ^ ( 1 - ^ ) 3 4 + v)(-l-v + v) 8 + v) 38 Chapter 3. Element Types and Shape Functions From (3.16)-(3.18), the following equations can be derived a dv = [J?]" a . dv . a dx a La y (3.20) 1 J (3.21) ds = dxdy = det[J ']dvdv E where [J ] is the Jacobian matrix of the transformation given by Ei - a fit " dv av a dv a x dx . dv (3.22) W ([**] [y ]) Ei dv . The integrals U$ and T$ in (2.39) can now be evaluated in the local coordinates. n n The matrix form will be used. [U ] and [T ] are the matrices corresponding to Ei Ei and T ' , respectively. These matrices are 8x8 square matrices. Some of their elements, E n however, may not be used in thefinalmatrix assembly because subscript m in U^ and n T ^ only refers to the unknown nodes in the element. [U ] and [T ] are given by E Ei Ei T a ^ = V~L V[/3] • V[/3] ^ r UEi JS Bi = — f 1 fi Ei [T ] Ei = in which J-i f dx = / f*Ei JS I Bi a 8y J I Ei T Ei Ei E E a [D <] = dx Ifl = [Jf]" a Ia J E y a Sy [/3] ds J (3.24) <T J^JP] [P}det[J <}dvdv T \L a dx (3.23) [D ) [D }det[J ]dvdv 1 J-i ( 1 . a dv [/?] a 8v . (3.25) As the integrals in (3.23) and (3.24) are related to the element node coordinates [x ] Ei and [y ], [U ] and [T ] cannot be evaluated once and for all. Also it would be too Ei Ei Ei complicated to integrate (3.23) and (3.24) analytically. Instead, numerical integration 39 Chapter 3. Element Types and Shape Functions formulas are used. With Gaussian quadrature formulas, [U ] and [T ] become Ei 1 [U ] Ei = Ei Ns Ns — EE^[^*(*i.'*)] [^*(«i^)]M^ («i^)] r , (3-26) VEi j=l =l k N N s [T ] Ei s = ^EE^^^>^)] ^i,^)]^[Jf(^,^)] r (3-27) 3=1 fc=l where {v^v^) gives the sampling point location in local coordinates, Wj and W). are the weighting factors related to Vj and i/*, respectively, and Ns is the number of sampling points in one direction. Vj and v\. are the sampling point locations in the Gaussian quadrature formula. These locations and associated weighting factors are given in Tab.3.2[5]. If the integrand in (3.23) or in (3.24) is a polynomial of order 2Ns — 1 for one variable (y or v), the integral associated with that variable can be exactly calculated with Ns sampling points. Similar to (3.13) and (3.14), F Ei h and F * are calculated from [T ]. E Ei Table 3.2: Locations of sampling points and weighting factors for Gaussian quadrature N s i Vj or Vj 0 2 3 5 3.3.2 2 1 0 ±0.7745966692414834 4 Wj ±0.3399810435848563 ±0.8611363115940526 0 ±0.5384693101056831 ±0.9061798459386640 8 9 5 9 0.6521451548625461 0.3478548451374539 0.5688888888888889 0.4786286704993665 0.2369268850561891 Quadratic triangular isoparametric element For the quadratic triangular isoparametric element shown in Fig. 3.6, all the previous equations, (3.15)-(3.17) and (3.20)-(3.25), are applicable except that the number of nodes is six instead of eight. The shape functions are the same as those of the second order simplex element, as indicated by the similarity between Fig. 3.1(b) and Fig. 3.6(b). Chapter 3. Element Types and Shape Functions 40 (0.0) {05,0) (1,0) V (b) local coordinates (a) global coordinates Figure 3.6: Quadratic triangular isoparametric element Using the node numbering scheme shown in Fig. 3.6 and assuming Ci = v a n d C2 = , v the shape functions can be derived from (3.2) and (3.4) as Pi = aoo = Ci(2& - 1) = v(2v - 1) 2 02 = "020 = C (2C - 1) = v(2v - 1) 2 03 = «oo2 = 2 Ca(2C - 1) = (1 - v 3 -2v- (3.28) 2v) 04 = "no = 4CiC = 4w 2 0s = "on = 4C2C3 = 4i/(l - v - v) 06 = aioi = 4CiC3 = 4v(l - v - u) The numerical integration formula is applied again to calculate [U ] and [T ] Ei Ei \U \ = ^E^il^(«i^i)] P*(«i»^)]^[^(«i.^M (3-29) [T ] = ^E^^i^i)] ^i^i)]^[^K^i)] (3-30) Ei Ei T r The locations and associated weighting factors are given in Tab.3.3 [11]. The error order term o(hk) indicates that the integral can be exactly calculated if the integrand is a 41 Chapter 3. Element Types and Shape Functions polynomial of order k — 1. Compared with simplex elements, fewer elements are needed if the above isoparametric elements are used to mesh regions with curved boundaries. Consequently, the number of nodes can be reduced. The procedure for evaluating the integrals of these isoparametric elements, however, is much more complicated. The integrals are solved numerically and have to be calculated within each element. Therefore, the node number reduction with isoparametric elements may not result in a reduction of computation time, as more time may be needed by the integral evaluations. F Ei k and F * are calculated from [T ]. E Ei h Table 3.3: Locations of sampling points and weighting factors for quadratic triangular element Error N Wi s i V15> 1 5 / 1 2 1 6 1 6 1 6 27 "96 25 96 25 96 25 96 V3> 3 / (0.47014206,0.05961587) (0.47014206,0.47014206) (0.05961587,0.47014206) (0.10128651,0.10128651) (0.79742699,0.10128651) (0.10128651,0.79742699) 0.066197075 0.066197075 0.066197075 0.06296959 0.06296959 0.06296959 ( 3 ' 3) (0,0.5) 3 (0.5,0) (0.5,0.5) V3> 3 / 4 (11 J.) M 5 > 15 / (2. 11) M5> 1 5 / (± 7 ±) 0(^2) o(h ) 3 o(h ) 4 0.1125 o(h ) 6 42 Chapter 3. Element Types and Shape Functions 3.4 Calculation of Integrals in the Loss-Energy Method In the loss-energy method for [Z] calculations, integrations (2.54) and (2.56) are used to calculate the power losses and stored magnetic energy. With the shape functions given in the preceding two sections, these integrations can be found. For high order simplex elements, in (2.54) becomes = £ j Jj^Jfc^ £ = ^ • nT,][J*J] (3-31) in which I refers to elements in conductor k, [J^] and [J^y)] are node value vectors r of current density distributions «/(«) and Jfa in element Ei, respectively. is the conjugate of J(kj)- Similarly, IOM,-,- in (2.56) is given by w Mij = E Re k=l / \ l = E Re fe SE,[A^] [T )[J^]\ Awfads) I S / JSB (3.32) T k=l \ I J in which [A^] is the node value vector of magnetic vector potential distribution A(ki) in element Ei. For quadratic isoparametric elements, = E E A f t ] fe=l Mij w = and WM^ are respectively given by l r ^ ' ] f t l (3-33) Ei a ERefE^-[4o] ^'][« S)]) r k=i \ i E 7 t ) (3-34) where [T ] is given either by (3.27) for quadrilateral elements or by (3.30) for triangular Ei elements. 3.5 General Procedures for [Z] calculations with F E M With the shape functions discussed in the preceding sections, all the entries in the final matrix given by (2.35) can be calculated. This section concentrates on the procedures for solving the final equations (2.35). 43 Chapter 3. Element Types and Shape Functions As mentioned in Section 2.4.2 that boundary nodes are numbered after unknown nodes, vector [A] can be partitioned as [Av] [A] = (3.35) [As] where [Av] = [A , A ,..., x 2 A] T N and [A ] = [A ,A ,, ,, AN ] . T B N+1 T N+2 Correspondingly, [U] + jw[T] is partitioned as [U]+MT} = (3.36) [D] [D ] B where [D] is an N x N sparse complex symmetric matrix, [D ] is an N x N complex B B matrix. With (3.35) and (3.36), (2.35) becomes ' [D] ~[F}' [hi] ' [Av] ' [Sc] . (3.37) . m+ [/«]. in which [F ] B = [IEI] = [ISI] = -MG ][F] T C (3.38) -[DB][A ] B MG ][F ] [A ] T C B B [IEI] and [7^2] are equivalent current vectors due to Dirichlet boundary conditions. In most cases of [Z] calculations, the boundary value vector [A ] is a zero vector. B In order to make use of the sparsity of [D], the node numbers are generally renumbered with certain algorithms [15] such that the matrix has the structure shown in Fig. 3.7. The shaded areas represent non-zero elements. [D] can now be factorized by algorithms dealing with banded symmetric complex matrices into [D] = [D ][Dv] L (3.39) Chapter 3. Element Types and Shape Functions 44 Figure 3.7: Matrix structure of the final equations for [Z] calculations [Di] is a banded lower triangular matrix and [Du] is a banded upper triangular matrix. Using the factors stored in [Di], (3.37) becomes '[Du] [F>] '' [Au] ' . [FH] [So] m . M (3.40) . [I] + [/«] . . in which [DL][F'\ = —[F] [DL}[I' ] = [IEI] EI (3.41) Because [Du] is an upper triangular matrix, it can be used directly to delete [FH] in (3.40), and the following equation can be derived '[Du] [F>] '' [Au] ' . PI [S'c]. . [ s) . WEI] (3.42) . [I] + [Ik*] . J in which [S'c] [Sc]-[F' ][F'] [I' 2\ [lE2]-[F' ][I' ] E [F' ][Du] H H H [FH] E1 (3.43) Chapter 3. Element Types and Shape Functions 45 [7] in (3.42) is left intact. It is the main excitation factor in the [Z] calculation because [/EI] and [IEI] are generally zero vectors. According to the discussions in Section 2.4 and Section 2.5, for an iV-conductor system (2.35) has to be solved N times in order to calculate the whole [Z], with only one conductor carrying current at a time. (3.42) can be used repeatedly to solve for [Js] and [Au] for different assignments to [I]. Once [I] is assigned, [Js] and [Au] are solved by [S' ][Js] c [I] + [iiJ (3.44) [Du][Au] The general procedure for [Z] calculations with FEM is summerized in the program flowchart shown in Fig. 3.8. 3.6 [Z] Calculation of an SC Coaxial Cable with F E M In this section, FEM is applied to calculate [Z] of a two-conductor SC coaxial cable. The numerical results are compared with analytical results given by (2.68) in Section 2.6. Optimum division and computation efficiency are studied and comparisons are made among different kinds of elements and different orders of simplex elements. Both the Js method and the loss-energy method discussed in Section 2.5 are used. The results show that for [Z] calculations of SC coaxial cables isoparametric elements are more efficient with respect to CPU time and storage requirements compared with simplex elements under the same accuracy. The loss-energy method gives the same results (up to eight digits) as the Js method. When solving the final symmetric banded complex linear equations, a partial pivotal element selection algorithm gives the same solutions as a variable bandwidth Choleski algorithm. Chapter 3. Element Types and Shape Functions 46 (Start) / Read data / Process for simplex elements I Renumber nodes 1 Form [U], [T], IF], [F ], and [G ] B c Assign / Form [ [[D][D ]] = [U]+j(o[T] andget B [ D] IFHI Aufl _ I" [ / ] "1 £1 [5 ] c J |[/,]J" [t/M/ ]J £ 2 Factorize [D] = [D u][D ] L I Assign [/ ] Solve [/ ] and U J s Calculate [Z] /"Plot fields and print [Z ] / Yes (Stop) Figure 3.8: Programflowchart for [Z] calculations with FEM 47 Chapter 3. Element Types and Shape Functions 3.6.1 Optimum divisions for SC coaxial cables Fig. 3.9(a) shows a two-conductor SC coaxial cable as an example for the [Z] calculation, and (b) shows its FEM solution region. Thefirststep in calculating [Z] with FEM is (a) cable data (b) solution region Figure 3.9: Geometry of a SC coaxial cable and its FEM solution region to mesh the solution region into elements. Thefinenessof the mesh will affect not only the accuracy but also the amount of computation. In general, afinemesh improves the accuracy but has more elements and nodes; therefore, it takes more CPU time and more storage. An optimum division scheme should give a mesh with the least possible number of nodes while maintaining a certain accuracy. For SC coaxial cables there are two relevant factors for thefinenessof the mesh: divisions along the radial direction and the span angle 0 shown in Fig. 3.9(b) Radial direction division Because of axial symmetry, only a wedge-shaped region as shown in Fig. 3.9(b) is needed for the [Z] calculation. A small suitable 6 can be used when studying the division scheme along the radial direction. The basic idea for meshing a solution region is to have fine elements at locations Chapter 3. Element Types and Shape Functions 48 where the field changes fast and to have large elements where the field changes smoothly. At high frequencies, the conductor currents will be concentrated in narrow regions near the conductor surfaces because of skin and proximity effects, and the magnitude of the current density decays quickly towards the centre of the conductors. In the region near the conductor surfaces, the field changes fast, and fine elements should be applied there in order to achieve accurate results. The depth of these regions beneath the conductor surfaces can be estimated by using the penetration depth (3.45) The divisions will then depend on 8, which is a function of frequency, conductor conductivity, and conductor permeability. The wedge-shaped region in Fig. 3.9(b) is enlarged in Fig. 3.10. The dashed lines in the figure are possible radial divisions. <£,j is the width of the jth division from a surface Figure 3.10: Radial divisions for SC coaxial cables of the ith conductor, dij < dij+i. For the inner most conductor, divisions proceed from the outer surface only, even if the conductor is hollow. For other conductors divisions proceed from both the inner and outer surfaces as shown in Fig. 3.10. The division width d^ is related to the penetration depth 8. This relationship can be 49 Chapter 3. Element Types and Shape Functions written as da = f Si (3.46) di where fd, is a division factor for the jth division and independent of Si. fdj is to be determined from numerical tests. The criteria are the error in [Z] and the error in the J distribution. These two errors are related to each other. For the cable shown in Fig. 3.9, the division factors found are listed in Tab.3.4, where Table 3.4: Division factors / j . for the SC coaxial cable element siml 0.3 iso or sim2 1.15 sim3 2.25 sim4 3.1 sim5 4.1 sim6 5 Si. 0.5 2.8 — — — — 0.9 /a. 1.8 — — — — — — — — — — "iso" stands for quadratic isoparametric elements including quadrilateral and triangular elements, "siml" stands for the 1st order simplex element, "sim2" stands for the 2nd order simplex element, and so on. Based on the division factors in Tab.3.4 and the data given in Fig. 3.9(a), the division radii for isoparametric and for the 2nd order simplex elements are calculated and listed in Tab.3.5. Division radii are the radii of division lines. The original material boundaries will also be division lines. Table 3.5: Division radii for iso and sim2 /(Hz) 6 60 600 6000 60000 600000 0 0 0 0 0 0 12 12 8.87 8.601 10.925 11.66 18 18 12 11.01 11.687 11.901 22 22 18 12 12 12 division radii (mm) 24 24 22 24 22 18 20 18 19.079 20.921 18 18.341 19.171 24 22 20.829 24 21.659 22 24 The meshes generated according to Tab.3.4 are used to calculate [Z], 0 = 6° is used in the calculation. The [R] and [L] values are listed in Tab.3.6. "ana" stands for Chapter 3. Element Types and Shape Functions 50 Table 3.6: [R] and [L] of a two-conductot SC coaxial cable R (fi/km) /(Hz) 6 60 600 6000 60000 600000 element Rn ana iso siml sim2 sim3 sim4 sim5 sim6 ana iso siml sim2 sim3 sim4 sim5 sim6 ana iso siml sim2 sim3 sim4 sim5 sim6 ana iso siml sim2 sim3 sim4 sim5 sim6 ana iso siml sim2 sim3 sim4 sim5 sim6 ana iso siml sim2 sim3 sim4 sim5 sim6 0.0388114 0.0388114 0.0388839 0.0388823 0.0388823 0.0388823 0.0388823 0.0388823 0.0417002 0.0417370 0.0418227 0.0418023 0.0417662 0.0417665 0.0417665 0.0417665 0.100575 0.100431 0.100645 0.100512 0.100645 0.100658 0.100643 0.100666 0.683376 0.682959 0.683443 0.682408 0.684061 0.682981 0.682998 0.683092 4.55387 4.54685 4.56176 4.55291 4.57740 4.55865 4.55860 4.55913 13.9902 13.9102 13.8688 13.9682 14.0963 14.0001 14.0028 14.0079 R\2 0.216122x10" B 0.215176x10"• 6 0.215349x10"• 6 0.214880x10"6 0.215736x10"6 0.215739x10" 6 0.215740x10--6 0.215740x10"• 6 0.216120x10"-4 0.215174x10"•4 0.215352x10"-4 0.214878x10"-4 0.215733x10"-4 0.215737x10"-4 0.215737x10"•4 0.215737x10"-4 0.00215824 0.00214979 0.00215072 0.00214678 0.00215449 0.00215453 0.00215453 0.00215453 0.190695 0.191040 0.191763 0.190726 0.190419 0.190437 0.190439 0.190439 1.70847 1.70634 1.71946 1.70905 1.71047 1.71048 1.71067 1.71055 5.11638 5.08706 5.08011 5.10549 5.12703 5.12097 5.12067 5.12252 L (jiH/km) R22 0.414466 0.414466 0.415225 0.415225 0.415225 0.415225 0.415225 0.415225 0.414477 0.414477 0.415235 0.415235 0.415235 0.415235 0.415235 0.415235 0.415564 0.415553 0.416317 0.416310 0.416320 0.416320 0.416320 0.416320 0.512251 0.512392 0.513133 0.513005 0.512891 0.512876 0.512877 0.512877 1.64196 1.64644 1.63998 1.64730 1.64424 1.64354 1.64278 1.64367 5.11640 5.08729 5.07796 5.10596 5.13941 5.12041 5.12070 5.12254 in 188.610 188.590 187.280 188.592 188.607 188.607 188.607 188.607 186.786 186.965 185.586 186.973 186.794 186.790 186.790 186.790 160.987 160.933 160.109 160.938 161.067 161.005 161.003 161.005 141.923 141.926 141.514 141.928 141.999 141.945 141.939 141.939 110.103 110.053 109.402 110.056 110.144 110.105 110.109 110.109 102.208 102.188 101.775 102.169 102.214 102.203 102.204 102.204 L22 36.1306 36.1314 36.1247 36.1261 36.1277 36.1282 36.1283 36.1283 36.1304 36.1312 36.1245 36.1259 36.1275 36.1280 36.1281 36.1281 36.1098 36.1141 36.1049 36.1085 36.1070 36.1074 36.1075 36.1076 34.2940 34.3098 34.3483 34.3066 34.2956 34.2972 34.2974 34.2975 21.5995 21.5830 21.4431 21.5819 21.6141 21.6007 21.5999 21.6004 18.7503 18.7471 18.7058 18.7377 18.7554 18.7487 18.7476 18.7475 29.4773 29.4758 29.2811 29.4689 29.4742 29.4747 29.4748 29.4749 29.4772 29.4758 29.2810 29.4688 29.4741 29.4746 29.4747 29.4748 29.4672 29.4676 29.2716 29.4605 29.4642 29.4647 29.4648 29.4648 28.5883 28.5986 28.4171 28.5928 28.5911 28.5885 28.5887 28.5887 21.6613 21.6629 21.5415 21.6573 21.6731 21.6620 21.6622 21.6619 18.7503 18.7467 18.7052 18.7375 18.7562 18.7487 18.7476 18.7475 Chapter 3. Element Types and Shape Functions 51 analytical results given by (2.68). Compared with analytical results the overall errors of the numerical results are less than 1%. The division radii listed in Tab.3.5 are also used for other frequencies with similar magnitudes. The overall errors of the numerical results in the frequency range from 0 to 1MHz are less than 1.2% for isoparametric elements, less than 1.8% for the 1st order simplex elements, and less than 1% for other order simplex elements. The corresponding current density distributions in the radial direction at 6kHz and 60kHz are plotted in Fig. 3.11 and Fig. 3.12, respectively. Span angle In practical problems, axial symmetry does not always exist to give the simple wedgeshaped solution regions shown in Fig. 3.9(b). Instead, a whole circular region as in Fig. 3.9(a) may have to be meshed. A larger span angle 0 results in a smaller number of elements and nodes. For straight-sided simplex elements, a larger 0 means a larger "misfit" of the element sides into the conductor shapes. This misfit will introduce errors. Curve-sided isoparametric elements can be fitted into the conductor shapes nicely, even with 0 larger than 90°. This enables isoparametric elements to maintain almost the same accuracy with large 0 as that with small 0. Fig. 3.13 gives the maximum errors in numerical results of [R] and [L] as functions of 0. The divisions along the radial direction are still based on Tab.3.4. For the 4th, 5th, and 6th order simplex elements the results are very similar to the 3rd order simplex element. From these results it can be seen that high orders in simplex elements do not really improve the accuracy when 0 is large. For isoparametric elements the maximum errors are less than 4% for [R] and 1.1% for [L] at 0 = 120°. 52 Chapter 3. Element Types and Shape Functions 0.020 0.020 • = ana V 0.015 = iso siml sim2 + = sim3 X sim4 A — d o.oio-i I V A o.oio- 0.005- o.ooo- •= 0.000\ [nr-JW" ~* 8 12 (mm) sim2 sim4 MH •'« M"igrJ IT -0.005 r (mm) 0.020 0.020 0.015 0.015 0.010 I — ana — —iso siml + = sim3 X 0.005 -0.005 • 0.015 • 7 -0 A — ana — —iso siml = sim2 + = aim 3 X = sim4 0.005 -0.005 -c -0.005 10.0 jag"* MP' <^ - 5.0- a 0.0 a = ana ^ = iso -» = siml A = sim2 + = sim3 x = sim4 -5.0 -10.0 -15.0 18 19 i 20 21 22 r (mm) Figure 3.11: Current density distribution in the SC coaxial cable at 6kHz Chapter 3. Element Types and Shape Functions Figure 3.12: Current density distribution in the SC coaxial cable at 60kHz 53 54 Chapter 3. Element Types and Shape Functions 46 0 80 30 (degree) 0 sim3 -o = 6 Hz = 60 Hz -•» = 600 Hz = 6 kHz -+ = 60 kHz -« = 600 kHz d ao 20 80 £2. 46 (degree) 16 a sim3 —a = 6 Hz —v = 60 Hz -~» = 600 Hz - -A = 6 kHz •-+ = 60 kHz —« = 600 kHz E 20 H o 9 46 0 (degree) 0 (degree) Figure 3.13: Maximum errors in [R] and [L] at different 0 Chapter 3. Element Types and Shape Functions 3.6.2 55 Computation efficiency: C P U time, storage, and pivoting From the preceding section it is seen that isoparametric elements can have larger 9 than simplex elements for the same accuracy. This will certainly reduce the number of elements and number of nodes in the mesh if a whole circular region needs to be meshed. Consequently, CPU time and memory storage will be reduced. There are other factors affecting the computation efficiency, such as bandwidth (BW) of thefinalequations depicted in Fig. 3.7, type of algorithm for solving thefinalequations, and element types. In order to achieve an overall comparison, a whole circular region is used in the comparison study. When the error limit is given, different element types with different 9 are applied to mesh the circular region. Fig. 3.14 gives the meshes of isoparametric and the 2nd order simplex elements at 60kHz with error limit of 2% and 15%, respectively. isoparametric 2nd order simplex (a) meshes with 2% error isoparametric 2nd order simplex (b) meshes with 15% error Figure 3.14: FEM meshes at 60kHz for different error limits In isoparametric element meshes, circles in division radii are also drawn for the purpose of misfit observation. The misfit in the mesh in Fig. 3.14(a) is almost unnoticeable, while it is very large in Fig. 3.14(b). Tab.3.7 lists some of the main computation parameters for different elements and different order of simplex elements at 2% error limit. The final matrix is a complex matrix in double precision. M and N are the number of elements 56 Chapter 3. Element Types and Shape Functions and number of unknown nodes, respectively. Tab.3.8 lists the corresponding CPU time for the calculation. The calculation is done on a VAX 11/750 having a speed of .6 ~ .7 million instructions per second withfloatingpoint accelerator. As a variable bandwidth Table 3.7: Storage and other parameters for different elements element type iso siml sim2 sim3 sim4 e (degree) 90 22.5 18 18 22.5 M N 32 464 300 220 176 89 225 581 961 1377 BW 17 20 119 206 319 matrix dimension 1088x1 3666x1 28474x1 70636x1 151385x1 storage (bytes) 83900 227084 1372376 3838008 7387372 Max errors (%) R L 0.85 1.87 1.65 1.81 0.06 1.59 0.52 0.25 — — Table 3.8: CPU time requirements for different elements element type ISO siml sim2 sim3 matrix formation 2.4 s (26.0%) 2.3 s (12.5%) 4.6 s ( 2.9%) 11.0 s ( 1.9%) matrix factorization 1.8 s (19.5%) 6.7 s (36.8%) 130.4 s (82.1%) 506.7 s (88.9%) solution 0.5 s ( 5.4%) 1.6 s ( 8.5%) 10.4 s ( 6.5%) 26.7 s ( 4.7%) loss-energy method 3.2 s (34.3%) 3.5 s (18.9%) 7.1 s ( 4.5%) 11.9 s ( 2.1%) others 1.4 s (14.8%) 4.3 s (23.4%) 6.4 s ( 4.0%) 13.9 s ( 2.4%) total CPU time 9.3 s 18.4 s 158.9 s 570.2 s Choleski algorithm was used[15], the final matrix was stored as a one-dimensional array. For the 4th order simplex element the calculation was not completed and the CPU time is not available. Tab.3.8 shows that the loss-energy method for calculating [Z] from the field solution takes a substantial amount of CPU time. The results, however, are the same as those from the Js method (up to eight digits), which takes negligible time. The agreement between these two methods exists in all the [Z] calculation of the SC cable. The variable bandwidth Choleski algorithm is modified from the one applied to real matrices. It takes advantage of the symmetry and sparsity of the matrix. However, it does not select a pivotal element when factorizing the matrix. Partial pivoting was also 57 Chapter 3. Element Types and Shape Functions used in the solutions, and the results were the same as those without pivoting. The CPU time required by the pivoting algorithm is much higher, and the corresponding CPU time and storage requirements are listed in Tab.3.9. Table 3.9: Storage and CPU time requirements for pivoting element type iso siml sim2 sim3 sim4 matrix dimension 49x89 58x225 355x581 616x961 955x1377 storage (bytes) 136268 377228 4216872 12179448 26005772 CPU time factorization solution 0.7 s 2.8 s 10.7 s 2.2 s 403.9 s 28.9 s 103.1 s 1703.7 s — — For 15% error limit, similar comparisons are given in Tab.3.10 to Tab.3.12. A Pivoting algorithm again gives the same results as the variable bandwidth Choleski algorithm. In Table 3.10: Storage and other parameters for different elements element type ISO siml sim2 sim3 sim4 sim5 sim6 e (degree) 180 60 60 60 60 60 60 M N 16 174 90 66 66 54 54 45 85 175 289 517 661 955 BW 9 9 29 61 105 161 229 matrix dimension 303x1 611x1 3158x1 9628x1 27845x1 50056x1 99610x1 storage (bytes) 33928 55604 170384 454104 1230332 2214940 4335336 Max errors (%) R 10.10 9.75 14.06 14.30 14.14 13.70 13.43 L 2.09 8.06 6.45 6.14 5.81 5.59 5.35 this case the loss-energy method gives faulty results with isoparametric elements. The reason is that the two centre triangular isoparametric elements in Fig. 3.14(b) have 180° as their vertex angles. Based on the above comparisons isoparametric elements are more efficient than simplex elements with respect to CPU time and memory storage within the same error limit. Chapter 3. Element Types and Shape Functions 58 Table 3.11: CPU time requirements for different elements element type ISO siml sim2 sim3 sim4 sim5 sim6 matrix formation 1.2 s (28.2%) 1.0 s (18.3%) 1.1 s ( 8.9%) 2.1 s ( 4.5%) 4.8 s ( 2.6%) 8.6 s ( 2.0%) 17.2 s ( 1.5%) matrix factorization 0.5 s (10.6%) 0.9 s (16.7%) 6.2 s (48.2%) 32.7 s (71.9%) 154.0 s (83.2%) 367.5 s (87.1%) 1017.2 s (90.0%) solution 0.2 s ( 4.6%) 0.4 s ( 6.8%) 1.3 s (10.3%) 3.7 s ( 8.2%) 10.9 s ( 5.9%) 18.4 s ( 4.4%) 39.7 s ( 3.5%) loss-energy method 1.7 s (38.7%) 1.3 s (25.5%) 2.2 s (17.3%) 3.7 s ( 8.2%) 7.9 s ( 4.3%) 11.1 s ( 2.6%) 19.9 s ( 1.8%) others 0.8 s (17.8%) 1.7 s (32.7%) 2.0 s (15.2%) 3.3 s ( 7.2%) 7.6 s ( 4.1%) 16.1 s ( 3.8%) 36.8 s ( 3.3%) total CPU time 4.4 s 5.3 s 12.8 s 45.5 s 185.2 s 421.7 s 1130.8 s Table 3.12: Storage and CPU time requirements for pivoting element type iso siml sim2 sim3 sim4 sim5 sim6 matrix dimension 25x45 25x85 85x175 181x289 313x517 481x661 685x955 storage (bytes) 47080 79828 357856 1137000 3373948 6501100 13208376 CPU time factorization solution 0.6 s 0.3 s 1.1 s 0.5 s 11.8 s 2.2 s 66.9 s 6.6 s 333.0 s 23.1 s 897.6 s 44.1 s 2569.4 s 117.3 s 3.7 Summary In this chapter shape functions are discussed for high-order simplex elements and for quadratic isoparametric elements. The integral matrices [Q^] and [Ts] of the 5th and 6th order simplex elements are given in the exact integer form. This supplements the tables given in [7]. The general procedure for [Z] calculations with FEM is also discussed. [Z] of a two-conductor SC coaxial cable is calculated by FEM with both simplex elements and isoparametric elements. Division factors for meshing SC cables are obtained through numerical tests. The results show that isoparametric elements are more efficient in the [Z] calculation of SC coaxial cables than simplex elements. Under the same error Chapter 3. Element Types and Shape Functions 59 limit isoparametric elements can have larger span angles; consequently, the FEM mesh has fewer elements and nodes. For large span angles the accuracy cannot be improved by increasing the order of the simplex elements. The results also show that the lossenergy method discussed in Section 2.5.2 takes a large amount of CPU time and gives the same results (up to eight digits) as the Js method. In solving thefinalsymmetric banded complex equations the partial pivotal element selection algorithm gives the same solutions as the variable bandwidth Choleski algorithm. Chapter 4 Earth Region Reduction Technique for [Z] Calculation with F E M 4.1 Introduction Power cables may be buried directly in the earth or installed in ducts or tunnels underneath the earth surface. Being a conductor itself, the earth often serves as a return path for the unbalanced currents in the cable system. Undoubtedly, the earth has to be included in the parameter calculation of underground power cables. The earth is generally represented as a uniform half space imperfect conductor for the parameter calculation of underground cables. With this assumption a formula for the impedance of a shallowly buried SC coaxial cable can be derived from the field distribution of a buried current filament by applying the perturbation concept. It is uncertain, however, whether this formula can be applied to other systems with different structures, such as tunnel installed cable systems. It is also uncertain where the perturbation concept becomes invalid as the earth penetration depth becomes smaller and smaller. The above uncertainties can be studied by FEM. Because the earth region extends to infinity, it is impossible for FEM to solve the field in the whole region. Fortunately, the earth never appears as an independent conductor, and it always exists as a reference conductor in the cable system. Consequently, the field to be calculated in the earth is always created by a loop current between one of the conductors in the cable system and the earth, and is concentrated around the cable system. Therefore, a boundary can be set up in the earth for FEM at a location sufficiently far away from the cable, where the 60 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM61 field becomes negligible, and FEM can then be used to calculate the parameters of the cable system together with the earth. Such a boundary is referred to as afieldtruncation boundary, and FEM with such a boundary shall be called "conventional FEM" here. With the help of the penetration depth, thefieldtruncation boundary for FEM can easily be established. In the earth the penetration depth can be very large in the low frequency range due to the poor conductivity of the earth. As a result, the FEM solution region becomes very large. More elements are needed to mesh such a solution region, and the computation time will increase. This is a weak point for the conventional FEM with thefieldtruncation boundary. There are several techniques for FEM to handle problems with infinitely large regions. Among them are the ballooning technique and the singular element technique. The ballooning technique can handle problems with regions of known Green's functions [45], but may be difficult to apply to the earth impedance calculation. The singular element technique assumes an approximate function representing the decay pattern of the field and creates infinitely long element sides to cover the infinitely large region by using singular shape functions. It is difficult tofindsuitable decay functions for this technique, although it is possible in principle to apply it to cable systems. In this chapter, a technique is proposed to reduce the earth region when the earth penetration depth is large. It is based on the same perturbation concept used to derive the impedance formulas of directly buried SC coaxial cables from thefieldsolution of a currentfilament.The boundary established by this technique has non-zero values, while the conventional field truncation boundary has zero values. The non-zero boundary values and the earth return current surrounded by such a boundary are calculated from thefieldsolution of the equivalent currentfilament.The proposed technique creates a much smallerfixedsolution region, and computation time can be saved if the earth return current is calculated only once. It is much easier to mesh a small region than a very large Chapter 4. Earth Region Reduction Technique for [ZJ Calculation with FEM 62 one. For a deeply buried SC coaxial cable, the earth can be assumed to occupy the whole space, and thefieldin the earth can be solved theoretically. Therefore, such a cable is usedfirstto study the division schemes and the locations of thefieldtruncation boundary in the earth for the conventional FEM. The impedances of the cable are calculated by the FEM with thefieldtruncation boundary and with the new proposed technique. The results are compared with those of analytical formulas. The comparisons show that the conventional FEM gives accurate results if thefieldtruncation boundary is at a location three times the earth penetration depth away from the cable. For the new proposed technique, the accuracy of results depends on the ratios rj,/5 and r /S . The parameters e e e r-fc and r are the boundary radius and inner earth radius, respectively. For r =24mm, e e accurate results can be obtained if rt,/8 < 0.2. e The impedances of shallowly buried and tunnel installed SC coaxial cables are also calculated with the conventional FEM, as well as with the proposed technique. The numerical results show that r > 12S is required for the conventional FEM. Comparisons 0 e with the conventional FEM for shallowly buried cables show that Pollaczek's formula gives noticeable errors when the earth penetration is small. The results of a tunnel installed cable show that Pollaczek's formula can also be applied to such a cable by using an approximate r if the earth penetration is large. e 4.2 Analytical Formulas for [Z] of Buried SC Coaxial Cables For shallowly buried SC coaxial cables approximate formulas for the [Z] calculation can be derived from thefieldsolution of an equivalent currentfilamentby applying the perturbation concept. For deeply buried SC coaxial cables, the earth can be assumed to occupy the whole space, and the correspondingfieldsbecome axisymmetrical. This leads Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM to analytical formulas. Fig. 4.1(a) shows a shallowly buried SC coaxial cable. In order to calculate [Z] of the cable, the field distribution in the earth has to be found. The impedance formulas can be derived from the E value on the inner surface of the earth, on which point P sits as shown in Fig. 4.1(a). The earth always serves as a return path for the current unbalance in the cable, and the field inside the earth is caused by a loop current between the cable and the earth. This loop current is the sum of all the conductor currents in the cable. (a) real structure (b) equivalent current filament Figure 4.1: A shallowly buried SC coaxial cable No formula exists for finding the field distribution in the earth caused by the above loop current if the actual geometry of the cable is considered. To overcome this difficulty, an equivalent current filament shown in Fig. 4.1(b) is used to replace the original cable. The filament is located at the centre of the cable and is insulated from the earth. If the earth is assumed to be uniform in half space, with its surface being parallel to the buried current filament, the field distribution in the earth caused by the loop current between the filament and the earth can be solved analytically. Such field solutions were first given by Pollaczek [2] and later by Wedepohl et a/[10]. Applying the assumptions made in Section 2.2 to the structure in Fig. 4.1(b), with 63 Earth Region Reduction Technique for [Z] Calculation with FEM64 Chapter 4. x-y axes as shown in thefigure,the Efielddistributions are [2] [10] ^ jwu.J f°» -va-fc '«*+i/p2 _ j a ^ i f°°_e * os(sa)<fa e = C 7T E = e y>0 v 2TT 7o - K \ a 0 + ^ ^ a 2 + — - K 0 \PeJ 1 (4.1) /p2 — \pj + / Jo e± . a cos{xa)da + ^ 2 + i/ 2 a p y<0 (4.2) where = v/x + (y + /i) = y/ * + (y-hy £' 2 (4.3) 2 (4.4) x Pe = (4.5) E and J? are the Efieldsin the air and in the earth, respectively, and their detailed a e derivations are given in Appendix B. h is the burial depth of the cable. I is the current in thefilament.p is the complex penetration depth in the earth. Ko is the zero order e second kind modified Bessel function. fi is the permeability of the air. /x and p are, a e e respectively, the permeability and the resistivity of the earth. The penetration depth in a conductor is a measure of the field attenuation in the conductor. It can be used to approximately find the region within which the most significant part of thefieldexists. For convenience, real penetration depth S defined in (3.45) is often used which relates to p through S=y/2\p\. Because the earth is a poor conductor, its penetration depth will be quite large, especially in the low frequency range. With fi = / i and typical value of p = 10012m for the earth, the real penetration depth e 0 e in the earth S is 5033 m at 1 Hz and 5.03 m at 1 MHz. e Based on the perturbation concept, the impedance formulas for shallowly buried SC coaxial cables can be derived from (4.2). Considering the large differences between the Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM65 earth penetration depth and the dimensions of the cable structure, it is assumed that the field of the current filament in Fig. 4.1(b) would not be disturbed significantly if the filament were replaced by the original SC coaxial cable in Fig. 4.1(a). Therefore, E (xp, yp), the E value at point P on the inner surface of the earth as shown in Fig. 4.1(a), e can be calculated approximately from (4.2). The surface impedance at point P is defined as Z E (xp,y ) I e = E P w« ^ ^ + = + / Jo . + (yp + h) 2 j K u ^ y t e + (yp ~ h y y = cos(xp ct)da 1 / + (4.6) Strictly speaking, Z would be different at different locations on the inner surface of the E earth. The differences, however, are very small, and can be ignored. If xp = r and e yp = —h, (4.6) becomes , ( Z.= 2^r \ K o (- ) - K ( £ ± ^ \ \p J e 0 \ ) p e + f Jo ^ a + yja* + 1/pl cos(r«a)J J (4.7) r is the inner surface radius of the earth. Z is also called "earth return impedance." e E Equation (4.7) was first derived by Pollaczek, and shall be called "Pollaczek's formula." Once Z is known, [Z] of a buried SC coaxial cable can be calculated with (2.68), E except that ZEQ K in (2.72) has to be modified into ZEQ K = ZB K + ZD K + Z E (4.8) Z in (4.7) can be found by numerical integration. It should be noted from (2.69) and E (2.70) that every element in matrix (2.68) has ZEQ K as one of its components. Therefore, Z must be added to every element in [Z] of a SC coaxial cable. E Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM66 If a SC coaxial cable is buried deeply in the earth (h ^> S ), the earth can be assumed e to occupy the whole space, and the field becomes axisymmetrical again. For such a case the earth becomes a coaxial conductor with infinite outer radius, and its inner surface impedance will be derived from (2.65) by letting TB — • oo and = r . The resulting h e formula is Z e ~ 27rr e Mr^) ™> [ r is the inner surface radius of the earth. The above equation is to be used in (4.8) e together with (2.68) to calculate [Z] of deeply buried SC coaxial cables. There are some uncertainties when applying the above formulas to the [Z] calculation of practical underground cable systems. For example, for the tunnel installed SC coaxial cable shown in Fig. 1.1, (4.7) cannot be applied directly because r does not exist. Also, e it is not clear how the tunnel structure will affect the impedance calculation with respect to the insulation within the tunnel. Similarly, Z for shallowly buried SC coaxial cables e is derived by using an approximation with the perturbation concept, and it is uncertain when this concept is no longer applicable. These uncertainties can be studied numerically with FEM. Before FEM is applied to the [Z] calculation of a general underground cable system, it is first applied to the [Z] calculation of a deeply buried SC coaxial cable, where the field can be found analytically. This case is used to study the division schemes in the earth and the locations of the field truncation boundary for FEM. It is also used to develop a new technique for reducing the earth region in FEM solutions in Section 4.4. Chapter 4. 4.3 Earth Region Reduction Technique for [Z] Calculation with FEM 67 [Z] Calculation of Deeply Buried SC Coaxial Cables by Conventional F E M with a Field Truncation Boundary It is assumed that the two-conductor SC coaxial cable in Fig. 3.9(a) is buried deeply in the earth as shown in Fig. 4.2(a). To calculate [Z] of the cable numerically by conventional FEM, afieldtruncation boundary located at r& shown in Fig. 4.2(b) has to be established, which should enclose the most significant part of thefielddistribution. Because the earth is a very poor conductor, the division scheme for conductors discussed (a) a deeply buried cable (b) FEM solution region Figure 4.2: A deeply buried SC coaxial cable and its FEM solution region in Section 3.6 is not applicable to the earth, and new division schemes need to be found. Two questions to be answered in this section are therefore: where to locate the field truncation boundary and how to mesh the earth. For the system in Fig. 4.2(a) [Z] can be calculated analytically. Therefore, the numerical results from FEM can always be checked by comparing them with those from analytical formulas. As the earth serves as a return path or as a reference conductor, the general FEM solution procedure discussed in Section 3.5 needs to be modified slightly. Assuming a system has K+l conductors with the reference conductor being conductor K+l, [Z] will Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM be a K x K matrix. To calculate the jth column of [Z] , a loop current Ij between the jth conductor and the reference conductor is used to excite the field. The jth element and the (if+l)th element of [J] in (3.37) will be Ij and — respectively. Once the field is solved the elements in the jth column of [Z] are given by the following formula with the Js method Zij = T(JsJn where Js and Js { K+1 ductors, respectively. - Js K+1 /<TK+I) (i = 1,2,..., K) (4.10) are the source current densities in the ith and the reference con0% and <TK+I are the conductivities in the ith and the reference conductors, respectively. With the loss-energy method, [Z] is given by *** = ( E / c L *i = where { E R E Jjt ^ds]l{I ,I +I I ) R ( / C Ii (z,j = l , 2 , . . . , / 0 Ii W(V )}/( ^+W/>) 5 5 Ri 7 (*,j = l,2,...,2iQ (4.11) (4.12) and A(ki) are, respectively, the current density distribution and magnetic vector potential distribution caused by the loop current between the ith conductor and the reference conductor. For both methods the system has to be solved K times. To locate the optimum boundary is to find the minimum 77,. For the system in Fig. 4.2 the earth current enclosed by r can be calculated analytically. If a loop current of 1+j'O A e is assumed between the cable and the earth, with typical resistivity and permeability values as shown in Fig. 4.2(b), the earth return current I from the analytical formula is e shown in Fig. 4.3(a), as a function of rj, at 6kHz. S is 64.97m at this frequency. Similar e curves can be obtained for other frequencies. It can be seen that I eR wl A and I « 0 ei at 5S . This means that the most significant part of the field is enclosed by a boundary e at 5S for the deeply buried SC coaxial cable. e FEM is applied to calculate [Z] of the system in Fig. 4.2 for different r . Isoparametric 0 elements are used in the calculation, and the solution region is reduced to a wedged region 68 Chapter 4. 1 Earth Region Reduction Technique for [Z] Calculation with FEM 69 2 3 4 6 b/ e (a) enclosed earth current by r& at 6kHz r 6 1 2 3 4 9 b/ e (b) max errors in [R] and [L] r 6 Figure 4.3: Earth return current and maximum errors in [Z] for different r e similar to the one in Fig. 3.9(b) with 0 = 6° due to axial symmetry. Six frequencies are used: 6Hz, 60Hz, 600Hz, 6kHz, 60kHz, and 600kHz. To mesh the solution region, the division radii listed in Tab.3.5 are used for meshing the conductors, while the division pattern 10 , 10 ?, n n+ 10 5, n+ 10 n+1 discussed later in this section is used for meshing the earth. The maximum errors in R and L are plotted in Fig. 4.3(b). These maximum errors are selected from the errors in all the six elements of [R] or [L] at all the six frequencies. The figure shows that accurate results can be obtained for the deeply buried SC coaxial cable when r5=35 . e For the same system, p is varied from lOOOfhn to O.Olfim, and the results show that e the maximum errors in R and L always decrease to less than one percent when r{,=3£ . e Therefore, for [Z] calculations of deeply buried SC coaxial cables, 3£ can be used as the e location of a field truncation boundary for the conventional FEM. The optimum division in the earth means that the corresponding FEM mesh should achieve accurate results with the least possible number of elements and nodes. Due to the large difference between the size of the cable and the size of the earth region, it is impossible to have similar element sizes everywhere in the whole region. The general Earth Region Reduction Technique for [Z] Calculation with FEM 70 Chapter 4. practice is to use element sizes similar to those in the cable to mesh the areas in the earth region next to the cable. When meshing the areas in the earth region away from the cable, the element sizes will be increased gradually to the boundary at r . 0 Several division radius patterns are tried out for isoparametric elements with different earth resistivities. These patterns use several division radii in each decade from 10" to 10" . These patterns are listed in Tab.4.1. The first earth division radius is r and the +1 e last one is r =3S . Between r and r&, the earth is divided by division radii listed in 0 e e Tab.4.1. Table 4.1: Earth division radius patterns for isoparametric elements in the decade from 10" to 10" pattern division radii 1 10" 10" i 10" * 10" 10" 10"+2 10" 2 3 10" i x 10" \ x 10" 10 4 10" \ x 10" 10 5 10" 10" +1 + + +1 +1 +1 +1 +1 n+1 n+1 +1 The numerical results show that the first pattern in Tab.4.1 always achieves high accuracy. For the system in Fig. 4.2, the division radii in the earth given by this pattern at different frequencies are listed in Tab.4.2, and the resistances and inductances calculated Table 4.2: Earth division radii for isoparametric elements /(Hz) 6 60 600 6000 60000 600000 0.024 0.0464 0.1 46.4 100 215 0.024 0.0464 0.1 46.4 100 215 0.024 0.0464 0.1 46.4 100 215 0.024 0.0464 0.1 46.4 100 195 0.024 0.0464 0.1 46.4 61.6 0.024 0.0464 0.1 division radii (m) 0.215 0.464 1 464 1000 2150 0.215 0.464 1 464 1000 1949 0.215 0.464 1 464 616 0.215 0.464 1 2.15 4.64 4640 6164 2.15 4.64 10 21.5 2.15 4.64 10 21.5 2.15 4.64 10 21.5 10 21.5 0.215 0.464 1 2.15 4.64 10 21.5 0.215 0.464 1 2.15 4.64 10 19.5 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 71 from the corresponding meshes are given in Tab.4.3, together with the values calculated by the analytical formulas. The solution region has 9 = 6°. Table 4.3: [R] and [L] of the deeply buried SC coaxial cable L (mH/km) R (ft/km) /(Hz) 6 60 600 6000 60000 600000 element Ru ana iso ana iso ana iso ana iso ana iso ana iso 0.0447332 0.0447221 0.100918 0.100821 0.692750 0.691494 6.60507 6.59130 63.7665 63.6482 605.816 604.401 R\2 R22 0.0059220 0.420388 0.0059109 0.420377 0.0592392 0.473695 0.0591056 0.473561 0.594334 1.00774 0.593213 1.00662 6.11239 6.43395 6.09938 6.42073 60.9211 60.8546 60.8077 60.7478 596.942 596.942 595.578 595.578 L\\ 2.41400 2.41099 2.18191 2.17968 1.92586 1.92364 1.67654 1.67495 1.41446 1.41308 1.17633 1.17555 L12 L22 2.26152 2.25486 2.25853 2.25188 2.03126 2.02461 2.02884 2.02219 1.80098 1.79434 1.79882 1.79218 1.56891 1.56320 1.56733 1.56162 1.32596 1.32602 1.32461 1.32469 1.09287 1.09287 1.09210 1.09210 The first pattern is also tested with different earth resistivities and with different r , e and the results show that overall errors in the frequency range from 0 to 1MHz for all the elements in [R] and [L] are less than 1%. The resistivities used in the study are 100012m, 1012m, 111m, 0.112m, and 0.0111m, with the other parameters remaining the same. r is e varied among the values: 50mm, 100mm, 250mm, 500mm, and 1000mm with />=100flm, e while the other parameters remain the same. When r is larger than 50mm, additional e division radii are needed for the insulation between the outer conductor of the cable and the earth. The first pattern listed in Tab.4.1 emerges as a good choice in meshing the insulation. If a unit loop current is assumed for the system in Fig. 4.2 between the inner conductor of the cable and the earth, the corresponding J distributions in the earth calculated by FEM at different frequencies are shown in Fig. 4.4. The analytical J distributions in the earth are also plotted in the figure. Good agreements between the numerical solutions and the analytical solutions can be observed from the curves in the figure. The J distributions Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM72 r (m) r (m) Figure 4.4: J distributions in the earth at different frequencies Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 73 in the conductors of the cable at corresponding frequencies remain almost the same as those plotted in Fig. 3.11 and Fig. 3.12 because both cases have the same conductor currents and the same conductor division radii. The results also show that the errors in [Z] are directly proportional to the errors in J on the inner surface of the earth. By definition the earth return impedance for a deeply buried SC coaxial cable is obtained by dividing the inner surface E of the earth by the earth current. Therefore, errors in the inner surface J of the earth will be directly reflected in the earth return impedance. In meshing the solution region of the system in Fig. 4.2, different angles 9 are used for the isoparametric elements and for the second order simplex elements. Based on the results at six frequencies (6Hz, 60Hz, 600Hz, 6kHz, 60kHz, and 600kHz), the overall errors with isoparametric elements at 0 = 90° are approximately 1.2% for R and 0.2% for L. The same accuracy can only be achieved with the second order simplex element at 9 = 15°. The Js and the loss-energy methods for calculating [Z] from the field solutions always give the same results (identical to 8 digits) for all cases used in this section. It can be concluded from the results presented in this section that the impedances of deeply buried SC coaxial cables calculated by FEM with a field truncation boundary are sufficiently accurate provided that the boundary is at least Z5 away from the cables e and that the earth is meshed properly. This can serve as a guideline for choosing the field truncation boundaries of other types of underground cable systems where analytical solutions are not available for comparison purpose. Though the earth division pattern 10 , 1 0 3 , 1 0 § , 10 n n+ n+ n+1 gives accurate results, the number of elements in the earth is quite large, especially in the low frequency range. This increases computation time. Also, the large difference between the size of the cable and the size of the earth region make it difficult to mesh the whole solution region unless Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 74 a specific auto-mesh program is used. It is therefore worthwhile to search for ways of reducing the earth solution region. One such technique is proposed in the next section. 4.4 Earth Reduction Technique In Section 4.2 the formulas for a shallowly buried SC coaxial cable are derived with two approximations: the cable structure will not disturb thefielddistribution of an equivalent filament current located at the centre of the cable, and the field inside the inner surface of the earth is still axisymmetrical. The second approximation ignores the influence of the earth on thefieldinside the SC coaxial cable. It also ignores the differences among the E values at different points on the inner surface of the earth. These two approximations can be removed by applying FEM with field truncation boundary. Therefore, the influences of these approximations for shallowly buried SC coaxial cables can be studied. The method can also be used for underground cable systems with arbitrary structures, such as the tunnel installed SC coaxial cables shown in Fig. 1.1. From the discussions and results in the preceding section, it is clear that the conventional FEM requires an earth solution region which is very large compared with the dimension of the cables. This increases computation time and creates problems for general auto-mesh programs in meshing such a large region as well as the details around the cable. To overcome these problems, a technique based on the perturbation concept is proposed to reduce the earth solution region. If the earth penetration depth is much larger than the cable structure, the structure will only slightly disturb the field of a buried current filament located at the centre of the cable. Therefore, a small solution region can be used for FEM with non-zero boundary values and with partial earth return current enclosed by the boundary. The boundary Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 75 values and the partial earth return current are calculated approximately from the E field solution of the filament. The E field solution of a deeply buried current filament can be derived analytically, while that of a shallowly buried current filament has been derived by Pollaczek[2] and Wedepohl et a^lO]. In general, the boundary is set up at a distance which is large compared with the dimension of the cable structure. The resulting boundary will then only be located in the earth and the air. Assuming that Ef, and A\> stand for the E and A values on the reduced boundary, respectively, the following equations are used in the solutions A = - ~ H b : (4.13) in earth (4.14) 1 E Aj = — in air h Js ]U3(T e where Js and <r are the source current density and conductivity in the earth, respectively. e E\, is given by filament formulas and is proportional to the filament current. For a system of K+l conductors with conductor K+l as the reference conductor, if [EB] are the Eb values of boundary nodes in the FEM mesh, [AB] in (3.35) will become [A ] = JM B +M ^ ± i (4.15) where [ly] is a vector with NB elements. If the boundary nodes are numbered such that the first NB nodes are the boundary nodes in the earth, and the remaining NB—NB B B nodes are in the air, the first NB elements in [ly] will be unity and the rest of the B elements will be zero. [AB] becomes a partially known vector. [IEA and [IE ] 2 m (3.38) become fe] = -ID.] ( - J f f [/*] = MGc][F ] T B 1 + ( - ^ = [W - [ W « r „ (4-16) + M ^ ± l ) = [I ] + [S ]Js „ EI Cv K (4.17) Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM76 where [J*] = [D ][E ]/U») VBA = -[GC}[F ] [EB) B B T B [4.10) [F ] = [D ][1V]/UV<TK+I) [Sc ] = [Gc][F ] [lv]/<r i V v b T B K+ [Fv] and [Sc ] are vectors. Putting [J^,] and [IEJ] in (4.16) and (4.17) into (3.37) and v moving Js K+1 terms to the left side produces [D] -[FA] [FH] [Sc ] . ' " [Au] ' (4.19) . [J] + [IEA] . A [F] becomes [FA] if [Fv] is subtracted from its last column, and [Sc] becomes [Sc ] if A [Sc ] is subtracted from its last column. [I] always corresponds to a loop current. For Y a loop current Ij between conductor j and the earth, the jth and (iif+l)th elements in [I] will be Ij and J , respectively. J is the partial earth return current enclosed by the ep ep reduced boundary. Once (4.19) is solved for the loop current, (4.10) will be used to find [Z] . The loss-energy method is not applicable here. Although [AB] in (4.15) is only partially known due to unknown Js K+1 in the earth, the corresponding boundary nodes are still treated as havingfixedboundary values. They do not appear in thefinalequations of the solution as shown by (4.19). Their values will be evaluated by (4.15) after solving for Js K+1 For a deeply buried SC cable, the field solutions of both the original cable and the buried currentfilamenthave simple forms. Therefore, the proposed technique is first tested for such a cable. The cable shown in Fig. 4.2(a) is redrawn in Fig. 4.5(a), while the corresponding buried currentfilamentis drawn in Fig. 4.5(b). Ec(r) and Ep(r) in Fig. 4.5 are the Efieldsassociated with the cable and thefilament,respectively. They Chapter 4. t Earth Region Reduction Technique for [Z] Calculation with FEM77 f- - earth A / ^ S C cable / it V (a) original SC cable (b) current filament replacement Figure 4.5: Replacing a deeply buried SC cable with a current filament are given by E (r) c jW(l I K {r/p ) e = 0 e Mr) r e 2TT (r / )Ki(r /p ) = Pe -^K (r/ 2TT 0 e P e >r e (4.20) e r >0 ) (4.21) I is the loop current between the cable or the filament and the earth. The partial earth return currents between radii r and r for both cases in Fig. 4.5 are e ( IVe) \ T { 'IP*) r > r T W.) = M ^ ) - ^ ) = - (4.22) e t J ( r e / p e ) K i ( r e / p e ) ) ( i - ^ ^ (4.23) where I (r) = I (r)\ ^ F epc r = -1(1 - (r/p^Mr/p,)) Q (4.24) 7ir(r) is the earth return current enclosed by radius r with the filament. By defining a perturbation coefficient Cp as * = 7 K l Pe (?) <- > 4 25 \PeJ Ec(r) and Jep (r) can be respectively related to Ep(r) and Ieppi?) as c £7 (r) = c E^r)/^ r >r e (4.26) Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 78 I (r) = epc CpI (r) r>r epF (4.27) e The values of Cp at different |r«/p | are plotted in Fig. 4.6. From Fig. 4.6 it can be e seen that Cp is close to 1.0 when the earth penetration is large compared with r , since e l-r 0.989- ^ 0.998- 0.997- 0.9960.00 0.02 0.04 0.06 0.08 0.10 0.00 0.02 0.04 K/P.\ 0.06 0.08 0.10 \ ./p.\ r Figure 4.6: Perturbation coefficient Cp as a function of |r /p | e e the influence of the cable structure on the field distribution in the earth of the filament becomes negligible in this case. This confirms the basic assumption on which the proposed technique is based. For |r /p |=0.02, /=lMHz, /> =100£2m, and /z=/io> |p |=3.56m and e e e e e r =71.2mm. e By using the exact formulas in (4.20) and (4.22) to find [EB] in (4.18) and I for [I] in ep (4.19), [Z] should become reasonably accurate. The impedances of the cable in Fig. 4.2 are calculated with reduced earth regions with exact E\> and I for rfr=0.24m, lm, 2.5m, ep 10m, and 25m, and the results are listed in Tab.4.4. The results show that accurate [Z] is obtained in the low frequency range where r& S . When the frequency is high e enough such that S becomes comparable to r or even smaller than r , the final matrix e e 0 of the equations becomes ill-conditioned, and the results becomes erroneous. This may be due to the fact that [AB] is an unknown vector for the deeply buried cable although it is treated technically as a known vector. Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 79 Different earth resistivities and r are tested for the above case with the r& listed in t Tab.4.4, and the numerical results show that the errors are less than 1% in [R] for <S and in [L] for r& < 2S . p is varied among lOOOftm, lOflm, lflm, O.lfim, and O.Olfim e e with r =24mm, while r is varied among 50mm, 100mm, 250mm, 500mm, 1000mm, and e e 2000mm with />=100fim. e Table 4.4: [Z] of the deeply buried cable found from a reduced earth region R (n/km) / n (Hz) (m) Rn Rl2 #22 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 0.0447332 0.0447332 0.0447332 0.0447332 0.0447332 0.0447332 0.100918 0.100954 0.100955 0.100955 0.100955 0.100955 0.692750 0.692605 0.692607 0.692607 0.692606 0.692577 6.60507 6.60474 6.60466 6.60467 6.60393 6.60439 63.7665 63.7597 63.7593 63.7588 63.7098 63.8667 605.816 605.736 605.711 605.635 598.441 586.913 0.00592198 0.00592198 0.00592200 0.00592198 0.00592198 0.00592198 0.0592392 0.0592391 0.0592392 0.0592393 0.0592391 0.0592392 0.594334 0.594326 0.594326 0.594326 0.594325 0.594295 6.11239 6.11280 6.11274 6.11275 6.11204 6.11248 60.9211 60.9191 60.9189 60.9184 60.8736 61.0229 596.942 596.914 596.892 596.826 590.205 579.558 0.420388 0.420388 0.420388 0.420388 0.420388 0.420388 0.473695 0.473695 0.473695 0.473695 0.473695 0.473695 1.00774 1.00773 1.00773 1.00773 1.00773 1.00770 6.43395 6.43415 6.43409 6.43410 6.43340 6.43383 60.8546 60.8592 60.8590 60.8585 60.8137 60.9630 596.942 596.914 596.892 596.826 590.205 579.558 6 60 600 6000 60000 600000 L (mH/km) Xl2 Irn 2.41400 2.41323 2.41274 2.41247 2.41191 2.41164 2.18191 2.18138 2.18086 2.18059 2.18003 2.17976 1.92586 1.92508 1.92457 1.92430 1.92374 1.92347 1.67654 1.67583 1.67530 1.67503 1.67448 1.67421 1.41446 1.41369 1.41317 1.41290 1.41235 1.41131 1.17633 1.17559 1.17507 1.17480 1.17388 1.18241 2.26152 2.26079 2.26028 2.26001 2.25945 2.25918 2.03126 2.03055 2.03002 2.02975 2.02920 2.02892 1.80098 1.80027 1.79975 1.79948 1.79892 1.79865 1.56891 1.56821 1.56768 1.56741 1.56686 1.56659 1.32596 1.32523 1.32470 1.32443 1.32388 1.32287 1.09287 1.09215 1.09163 1.09136 1.09047 1.09834 2.25486 2.25413 2.25362 2.25336 2.25280 2.25253 2.02461 2.02390 2.02337 2.02310 2.02254 2.02227 1.79434 1.79362 1.79310 1.79283 1.79228 1.79200 1.56320 1.56250 1.56197 1.56170 1.56115 1.56088 1.32602 1.32531 1.32478 1.32451 1.32396 1.32295 1.09287 1.09215 1.09163 1.09136 1.09047 1.09834 € Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 80 The earth is meshed with the pattern 10 , 10 *, 10 *, 10 . The numerical results n n+ n+ n+1 also show that the impedance errors are sensitive to Eb and 1^ and are proportional to the errors in J on the inner surface of the earth. With the proposed technique, the exact E\, and 1^ for the reduced solution region are replaced by the approximate ones from the filament formulas in (4.21) and (4.23). Errors are introduced into E\, and 1^ by such a replacement as indicated by (4.26) and (4.27). It has been shown, however, that these errors become negligible if the earth penetration depth is much larger than the cable structure. The impedances of the cable in Fig. 4.2 calculated with reduced regions with approximate Eb and I ep from the filament field solutions for different r are given in Tab.4.5. 6 I is calculated by numerical integration based on Ep(r). The FEM solution meshes are ep used for the numerical integration as well. With r =24mm, the results show a similar pattern as with the exact Eb and I : e ep accurate impedances are obtained for the r listed in Tab.4.5 if 8 is much larger than e e r , and erroneous impedances are obtained if S is comparable with r . For the r listed 6 e b b in Tab.4.5, if the earth resistivity is varied among lOOOfim, lOfim, lfim, O.lfim, and O.Olfhn with the same r , the errors are less than 1% in [R] for r < 0.2S and in [L] for e 0 e r < 0.5£ . b e For the same S , the larger the r , the larger the ratio r /S , and the larger the errors e in Eb and I ep e e e given by the filament formulas. Therefore, for a large r , the ratio r /S e e e has to be less than a certain value in order to achieve accurate results with the proposed technique. By varying r among 50mm, 100mm, 250mm, 500mm, 1000mm, and 2000mm e with />=100ftm, the results show that the errors are less than 1% in [R] for r /S < 0.012, e e and in [L] for r /S < 0.055. e e The frequency range within which the proposed technique achieves high accuracy Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 81 is mainly determined by the relationship between r& and 8 for small r , and by the e e relationship between r and 8 for large r . For the cable with />=100fhn, the ratio r /8 e e e e e becomes more influential in determining the frequency range if r > 250mm. e Table 4.5: [Z] found with the proposed technique based on Ep f (Hz) 6 60 600 6000 60000 600000 R (ft/km) n (m) Rii ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 ana 0.24 1.00 2.50 10.00 25.00 0.0447332 0.0447332 0.0447332 0.0447332 0.0447332 0.0447332 0.100918 0.100954 0.100955 0.100955 0.100955 0.100955 0.692750 0.692605 0.692607 0.692607 0.692606 0.692575 6.60507 6.60480 6.60470 6.60472 6.60399 6.60443 63.7665 63.7646 63.7643 63.7633 63.7145 63.8656 605.816 606.086 606.059 605.980 598.766 586.921 R\2 0.00592198 0.00592197 0.00592198 0.00592198 0.00592198 0.00592198 0.0592392 0.0592391 0.0592391 0.0592391 0.0592391 0.0592391 0.594334 0.594326 0.594326 0.594326 0.594324 0.594292 6.11239 6.11285 6.11278 6.11280 6.11210 6.11252 60.9211 60.9237 60.9235 60.9226 60.8780 61.0218 596.942 597.237 597.214 597.146 590.504 579.565 L (mH/km) R22 0.420388 0.420388 0.420388 0.420388 0.420388 0.420388 0.473695 0.473695 0.473695 0.473695 0.473695 0.473695 1.00774 1.00773 1.00773 1.00773 1.00773 1.00770 6.43395 6.43420 6.43414 6.43415 6.43345 6.43387 60.8546 60.8638 60.8637 60.8627 60.8181 60.9617 596.942 597.237 597.214 597.146 590.504 579.566 L u 2.41400 2.41322 2.41274 2.41247 2.41191 2.41164 2.18191 2.18138 2.18086 2.18059 2.18003 2.17976 1.92586 1.92508 1.92457 1.92430 1.92374 1.92347 1.67654 1.67583 1.67530 1.67503 1.67448 1.67421 1.41446 1.41369 1.41317 1.41290 1.41235 1.41132 1.17633 1.17556 1.17504 1.17477 1.17386 1.18240 L12 L22 2.26152 2.26079 2.26028 2.26001 2.25945 2.25918 2.03126 2.03055 2.03002 2.02975 2.02920 2.02893 1.80098 1.80027 1.79975 1.79948 1.79892 1.79866 1.56891 1.56821 1.56768 1.56741 1.56686 1.56659 1.32596 1.32522 1.32470 1.32443 1.32388 1.32287 1.09287 1.09213 1.09160 1.09133 1.09045 1.09833 2.25486 2.25413 2.25362 2.25336 2.25280 2.25253 2.02461 2.02390 2.02337 2.02310 2.02254 2.02227 1.79434 1.79362 1.79310 1.79283 1.79228 1.79201 1.56320 1.56250 1.56197 1.56170 1.56115 1.56088 1.32602 1.32530 1.32478 1.32451 1.32396 1.32295 1.09287 1.09213 1.09160 1.09133 1.09045 1.09833 If r =500mm with p—lOOilm, the minimum 8 will be 41.7m, and the corresponding e e frequency is /=p /(7r/z i5 )«15kHz. Therefore, the frequency range for the proposed 2 e e Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM technique is from 1Hz to 15kHz, and rj, is assigned a fixed value, say 2.5m, for this frequency range. For the conventional FEM, r >36 is required for the deeply buried 0 e SC cable, and the corresponding rj, for the above frequency range will be from 15,100m to 125m. Thus the proposed technique achieves the goal of reducing the earth solution region in the low frequency range in this case. 4.5 [Z] Calculations of Shallowly Buried SC Coaxial Cables with F E M In this section the conventional FEM and the proposed technique discussed in the preceding section are applied to the [Z] calculation of shallowly buried SC coaxial cables. The impedances from Pollaczek's formula are compared against those from the FEM for different earth resistivities and for different r . The impedances of a SC cable with an e arbitrary structure are also calculated with FEM and with Pollaczek's formula. The numerical results show that for the conventional FEM a solution region with r > 0 12S will give reasonably accurate results. Based on the comparison with the conventional e FEM it is shown that accurate results can be obtained with the proposed technique for different p and r . Good agreements are observed among thefielddistributions around e e the cable from the conventional method, from the proposed technique, and from the approximate analytical solution. The comparisons between Pollaczek's formula and the FEM show that accurate results are obtained with Pollaczek's formula when r /S is small. e e For r =24mm, the maximum differences between the results from Pollaczek's formula and e those from the FEM are less than 1% in R if r /S < 0.03 and in L if rJ8 <0.095. For e e e r„=1000mm, the maximum differences are less than 1% in R if r /S <0.018 and in L if e e e/^ <0.154. The numerical results for the cable with an arbitrary structure show that r e Pollaczek's formula can still be used with an approximate r . e 82 Chapter 4. 4.5.1 Earth Region Reduction Technique for [Z] Calculation with FEM 83 Determination of the solution region for the conventional F E M in [Z] calculations of shallowly buried cables With a shallowly buried SC coaxial cable, there is no exact field solution which could be used as a comparison base. The solution region for the conventional FEM is, therefore, determined by an iterative procedure. The boundary radius will be increased gradu- ally, and the difference between results from two consecutive iterations should become smaller as rj becomes larger. This difference is used as a criterion in the solution region determination. When the earth penetration is large, Pollaczek's formula (4.7) could also be used as a reference in the solution region determination. It is assumed that the SC coaxial cable shown in Fig. 3.9 is buried 1.5m beneath the earth surface. Therefore,fe=1.5mand r =24mm in Fig. 4.1(a). It is also assumed that e />=100f2m and / z e e =/i . 0 Once 77, is chosen, the earth is divided with the isoparametric elements in the pattern 10 , 10 3, 10" t, 10 n n+ + n+1 discussed in Section 4.3. The mesh for rj,=2£ at 6kHz is shown in Fig. 4.7. e (a) entire mesh (b) detailed mesh around the cable Figure 4.7: FEM mesh at 6kHz for r =2S b e The impedances calculated with r& being equal to 55 , 106,., and 15S , and with e e Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM Pollaczek's formula are listed in Tab.4.6. It can be seen that for [L] there is a good agreement between the results from the conventional FEM and from Pollaczek's formula at rh=5S . For [R] good agreement exists at r =15S . e b e Table 4.6: [Z] of the shallowly buried SC coaxial cable from the conventional FEM / R (fi/km) n (Hz) 56, 106, 156, 6 Pollaczek 56, 106, 156, 60 Pollaczek 56, 106, 156, 600 Pollaczek 56, 106, 156, 6000 Pollaczek 56, 106, 156, 60000 Pollaczek 600000 56, 106, 156, Pollaczek Rll Rl2 0.0445242 0.0446705 0.0447053 0.0447405 0.0989151 0.1004902 0.1008098 0.1011468 0.677598 0.692109 0.695547 0.699820 6.59260 6.74617 6.77707 6.81484 67.3478 68.6650 68.9727 69.3549 698.632 709.490 711.493 714.613 0.00570694 0.00585330 0.00588803 0.00592928 0.0571899 0.0587650 0.0590845 0.0594682 0.579357 0.593867 0.597305 0.601403 6.10068 6.25425 6.28515 6.32215 64.5070 65.8243 66.1319 66.5095 689.808 700.667 702.670 705.740 L (mH/km) R-22 0.420238 0.420384 0.420419 0.420395 0.471710 0.473285 0.473605 0.473924 0.992824 1.007335 1.010773 1.014809 6.42211 6.57568 6.60658 6.64371 64.4471 65.7643 66.0720 66.4430 689.808 700.667 702.670 705.740 L\\ 2.50775 2.50897 2.50916 2.51380 2.27548 2.27677 2.27695 2.28130 2.01881 2.01996 2.02014 2.02392 1.76564 1.76670 1.76683 1.77047 1.49249 1.49300 1.49306 1.49589 1.22252 1.22201 1.22186 1.22453 L\2 2.35534 2.35656 2.35675 2.36132 2.12469 2.12598 2.12615 2.13064 1.89399 1.89514 1.89532 1.89904 1.65802 1.65908 1.65921 1.66284 1.40402 1.40453 1.40458 1.40739 1.13908 1.13856 1.13842 1.14107 l>22 2.34869 2.34990 2.35010 2.35467 2.11803 2.11932 2.11950 2.12399 1.88734 1.88849 1.88867 1.89240 1.65231 1.65337 1.65350 1.65713 1.40410 1.40461 1.40466 1.40745 1.13908 1.13856 1.13842 1.14107 The earth return current enclosed by the boundary is calculated from E in (4.2) with e numerical integration over the FEM mesh, and its real and imaginary parts are plotted in Fig. 4.8(a). A loop current of 1+j'O A is assumed between the centre conductor of the cable and the earth. Comparing with Fig. 4.3 for the deeply buried SC cable, the earth current approaches its final value of 1+j'O much slower as rj, increases. In Fig. 4.8(b) is the plot of the maximum differences. The maximum differences between the [Z] calculated with rf,=nS and the [Z] calculated with r(,=(n-l)5 is denned as the nth maximum e e difference. It is taken from all the elements in [Z] in the frequency range from 1Hz to 84 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM85 rb/ e (a) enclosed earth current by r at 6kHz b/ e (b) max differences in [R] and [L] s r 0 6 Figure 4.8: Earth return current and maximum differences in [Z] for different r 0 1MHz. Fig. 4.8 shows that the maximum difference in [R] decreases slowly. It reaches 0.07% at r =15S . The maximum difference in [X] reaches 0.027% even at r =6S . 0 e 0 e If the nth maximum difference is defined as the one between r6=(ra-3).5 and r =nS e 0 e and if 0.5% is chosen as the threshold value of the maximum difference in [R] for stopping the iteration, r =12S will be the iteration result because the maximum difference in [R] 0 e is 0.55% at 12S and .27% at 15S . The numerical results show that the maximum e e difference in [R] between 12^ and 158 maintains the value 0.27% at different p and e e r . For p =100flm and r =24mm, the differences between the results found with the e e e FEM having r =12S and those found with Pollaczek's formula are less than 1.1% in 0 e [R] and less than 0.25% in [L] in the frequency range from 1Hz to 1MHz. Therefore, for the conventional FEM, a solution region with rj, >12S will give reasonably accurate e impedances for the shallowly buried SC cables. Chapter 4. 4.5.2 Earth Region Reduction Technique for [Z] Calculation with FEM 86 Application of the proposed technique to [Z] calculations of shallowly buried SC cables The proposed technique discussed in Section 4.4 is applied to the same cable used in the preceding subsection with h=1.5m, r =24mm, />=100fhn, and p =po. Four different r e e e 0 are used: 2.5m, 5m, 10m, and 25m. The meshes are similar to those in Fig. 4.7. The impedances found with the proposed technique at different r are listed in Tab.4.7. 0 "conventional" in the table represents the results found with the conventional FEM, and "Pollaczek" represents the results given by Pollaczek's formula. Good agreement can be observed among the three approaches. The partial earth return current I ep used in the calculations at 60Hz, 6kHz, and 600kHz are listed in Tab.4.8. They correspond to a l+j'O A loop current between the cable and the earth. If the earth return current is l+j'O A at 6kHz, the corresponding Efielddistributions in the air and the earth, given by (4.1) and (4.2), respectively, are plotted in Fig. 4.9. It can be seen that E is very smooth, except that it becomes singular at the location of the equivalent current filament. This smoothness makes the numerical integration easier. For the same earth current and frequency, the corresponding J contour lines in the earth from the three approaches are plotted in Fig. 4.10. "proposed" represents the results found with the proposed technique with r =5m, and "analytical" represents the 0 results derived from (4.2). The plotting area is 2.5m by 3.5m. Good agreement in the J distributions among the three approaches can be observed from Fig. 4.10. The efficiency of the proposed technique is compared with that of the conventional FEM. For the above calculations, storage and CPU time requirements for the proposed technique with rj,=5m are listed in Tab.4.9 and Tab.4.10, respectively. The CPU time is taken from a VAX-750. I ep represents the time in calculating the earth return current with numerical integration. M, N, and BW are the same as in Subsection 3.6.2. Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM Table 4.7: [Z] of the shallowly buried SC coaxial cable from the proposed technique n (Hz) (m) R\\ conventional 2.5 5 10 25 Pollaczek conventional 2.5 5 10 25 Pollaczek conventional 2.5 5 10 25 Pollaczek conventional 2.5 5 10 25 Pollaczek conventional 2.5 5 10 25 Pollaczek conventional 2.5 5 10 25 Pollaczek 0.0447186 0.0447465 0.0447465 0.0447465 0.0447465 0.0447405 0.100950 0.101193 0.101193 0.101193 0.101193 0.101147 0.696866 0.699635 0.699635 0.699636 0.699640 0.699820 6.79056 6.81436 6.81438 6.81443 6.81482 6.81484 69.0895 69.3455 69.3470 69.3523 69.3739 69.3549 712.336 714.293 714.421 714.270 713.826 714.613 6 60 600 6000 60000 600000 L (mH/km) R (ft/km) / R\2 0.00590137 0.00592928 0.00592928 0.00592928 0.00592928 0.00592928 0.0592246 0.0594680 0.0594680 0.0594680 0.0594681 0.0594682 0.598624 0.601394 0.601394 0.601394 0.601399 0.601403 6.29865 6.32245 6.32246 6.32251 6.32290 6.32215 66.2488 66.5047 66.5063 66.5115 66.5332 66.5095 703.512 705.469 705.597 705.447 705.003 705.740 #22 In 0.420432 0.420460 0.420460 0.420460 0.420460 0.420395 0.473745 0.473988 0.473988 0.473988 0.473988 0.473924 1.01209 1.01486 1.01486 1.01486 1.01487 1.01481 6.62008 6.64388 6.64389 6.64394 6.64433 6.64371 66.1888 66.4448 66.4463 66.4516 66.4732 66.4430 703.512 705.469 705.597 705.447 705.003 705.740 2.50922 2.51230 2.51207 2.51172 2.51144 2.51380 2.27700 2.28000 2.27977 2.27943 2.27915 2.28130 2.02019 2.02243 2.02220 2.02185 2.02157 2.02392 1.76687 1.76903 1.76880 1.76846 1.76818 1.77047 1.49306 1.49441 1.49418 1.49383 1.49350 1.49589 1.22179 1.22308 1.22284 1.22235 1.22215 1.22453 2/12 L22 2.35681 2.35015 2.35989 2.35323 2.35966 2.35300 2.35931 2.35266 2.35903 2.35238 2.36132 2.35467 2.12620 2.11955 2.12921 2.12255 2.12898 2.12232 2.12863 2.12198 2.12835 2.12170 2.13064 2.12399 1.89537 1.88872 1.89761 1.89096 1.89738 1.89073 1.89703 1.89039 1.89675 1.89011 1.89904 1.89240 1.65925 1.65354 1.66142 1.65570 1.66119 1.65547 1.66084 1.65513 1.66056 1.65485 1.66284 1.65713 1.40459 1.40467 1.40594 1.40602 1.40571 1.40579 1.40536 1.40544 1.40503 1.40511 1.40739 1.40745 1.13835 1.13835 1.13964 1.13964 1.13940 1.13940 1.13891 1.13891 1.13871 1.13871 1.14107 1.14107 Table 4.8: I found with numerical integrations for the proposed technique ep Re(J ) and Im(/ ) (p.u.) / ep (Hz) 60 6000 Re(J „) Im(I„) Re(/ ) Im(I«„) e ep 600000 r(,=2.5m 0.583064xl00.486857xl00.598188 xlO" 0.314763xl0" 0.641595xl00.130791 5 4 3 2 x Im(J„) ep rj,=5m 0.233561 xlO" 0.174161 xlO" 0.241264xl0" 0.105046X100.233091 0.310736 4 3 2 1 rj,=10m 0.936563xl0" 0.603340xl00.973100xl0" 0.324944 x 10" 0.646749 0.434091 4 3 2 1 r =25m 0.589085xl00.302400xl00.598538X10" 0.125884 1.005062 0.895586X10tt 3 2 1 1 87 JR (A/km ) 2 (A/km ) 2 J/.=250.20+17.26-(t-l) conventional J =250.95+17.25 (t-l) Horizontal distance (m) A proposed (r =5m) 0 analytical Figure 4.10: J distributions in the earth at 6kHz from the three approaches Chapter 4. Earth Region Reduction Technique for [Z] Calculation vrith FEM Table 4.9: Storage and other parameters for the proposed technique (rj,=5m) / (Hz) 60 6k 600k M N 283 343 403 88 106 124 BW 36 36 36 matrix dimension 4351x1 5356x1 6361x1 storage (bytes) 302064 365976 429888 Table 4.10: CPU time requirements for the proposed technique (r =5m) 4 matrix formation 6.8 s ( 8.6%) 8.7 s ( 9.9%) 10.1 s (11.3%) / (Hz) 60 6k 600k matrix factorization 9.9 s (12.4%) 11.8 s (13.4%) 13.8 s (15.4%) solution 8.0 s (10.1%) 7.7 s ( 8.7%) 8.1 s ( 9.1%) others 4.2 s ( 5.3%) 5.1 s ( 5.7%) 5.8 s ( 6.5%) 50.4 s (63.5%) 55.1 s (62.3%) 51.8 s (57.8%) total total-/,,, 79.3 s 28.9 s 88.4 s 33.3 s 89.6 s 37.8 s For the conventional FEM with r =12S the storage and CPU time requirements for 0 e) the same problem are listed in Tab.4.11 and Tab.4.12, respectively. Table 4.11: Storage and other parameters for the conventional FEM (rt,=126 ) e f (Hz) 60 6k 600k M N 463 463 463 142 142 142 BW 32 32 33 matrix dimension 9196 x 1 9196x 1 9259x 1 storage (bytes) 507600 509328 512064 Tab.4.9-Tab.4.12 indicate that less storage is needed for the proposed technique than for the conventional FEM due to a smaller solution region. As the frequency becomes large, the storage requirement becomes similar for both approaches. The total CPU time with the proposed technique is much higher than with the conventional FEM. On average, about 60% of the total time is spent on the 1^ calculation. As long as the earth parameters p , p and r , as well as the locations of the equivalent current filaments do e t% e not change, I only needs to be calculated once. By deducting the time for calculating ep I from the total CPU time, the proposed technique will require less CPU time than the ep conventional FEM. 89 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 90 Table 4.12: CPU time requirements for the conventional FEM (r =12S ) 0 f (Hz) 60 6k 600k matrix formation 11.3 s (25.8%) 11.6 s (27.0%) 11.9 s (27.2%) matrix factorization 21.9 s (50.0%) 20.7 s (48.4%) 21.3 s (48.8%) solution 4.1 s ( 9.4%) 4.0 s ( 9.3%) 3.9 s ( 9.0%) others 6.5 s (14.7%) 6.5 s (15.2%) 6.5 s (15.0%) e total 43.8 s 42.8 s 43.6 s If p is varied with r remaining at 24mm, the maximum differences between the ree e sults from the conventional FEM and from the proposed technique are listed in Tab.4.13. They are chosen from all the elements in [Z] in the frequency range from 1Hz to 1MHz. Table 4.13: Maximum differences in [Z] with the proposed technique at different p e Pe (fim) 1000 100 10 1 0.1 rj=2.5m in R in L 1.05% 0.15% 1.05% 0.15% 1.05% 0.14% 1.02% 0.14% 1.02% 0.13% rj=25m r =10m r»=5m in R in L in R in L in R in L 1.05% 0.14% 1.05% 0.13% 1.05% 0.12% 1.05% 0.14% 1.05% 0.12% 1.05% 0.11% 1.05% 0.13% 1.05% 0.12% 1.05% 0.10% 1.02% 0.13% 1.02% 0.11% 1.02% 0.09% 1.02% 0.12% 1.02% 0.10% 1.02% 0.08% 6 If r is varied with p remaining at lOOflm, the corresponding maximum differences are e e listed in Tab.4.14 Table 4.14: Maximum differences in [Z] with the proposed technique at different r r, (mm) 50 100 250 500 1000 rj=2.5m rj,=5m in R in L in R in L 1.05% 0.15% 1.05% 0.14% 1.05% 0.15% 1.05% 0.14% 1.52% 0.15% 1.23% 0.14% 4.11% 0.62% 2.89% 0.53% — — 7.99% 1.72% r =10m in R in L 1.05% 0.12% 1.05% 0.12% 1.13% 0.13% 1.91% 0.45% 4.56% 1.42% 6 e r»=25m in R in L 1.05% 0.11% 1.05% 0.11% 1.05% 0.11% 1.23% 0.15% 2.20% 0.48% From the results presented in this subsection it can be seen that for small r /8 e e accurate impedances of shallowly buried SC cables can be obtained with the proposed technique. With this technique, the solution region in the earth becomes small and fixed. In contrast, the solution region with the conventional FEM varies with the frequency in Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 91 order to save time, and becomes very large at low frequencies. Small regions are easier to mesh than larger ones. The proposed technique takes more CPU time compared to the conventional FEM, but if I 4.5.3 is calculated only once, it will require less time. ep Comparisons between analytical results and F E M results for shallowly buried SC cables The conventional FEM is used to investigate the errors of Pollaczek's formula under different conditions. With the same assumption that the cable in Fig. 3.9(a) is buried 1.5m beneath the earth, the impedances of the cable are calculated with both the conventional FEM and Pollaczek's formula for different p and r . The maximum differences between e e the impedances from the two approaches are listed in Tab.4.15, within the frequency range of 1Hz to 1MHz. Table 4.15: Maximum differences in [Z] with Pollaczek's formula at different p Pe (Clm) 1000 100 10 1 0.1 r,=24mm in R in L 0.48% 0.22% 0.48% 0.24% 0.48% 0.25% 2.56% 0.24% 9.66% 3.29% r«=100mm r =250m in R in L in R in L 0.48% 0.22% 0.64% 0.24% 0.97% 0.22% 2.80% 0.23% 3.28% 0.25% 8.85% 1.57% 12.36% 2.91% 20.29% 8.22% 20.34% 15.98% 47.96% 15.47% e e r =1000mm in R in L 4.35% 0.26% 13.35% 2.20% 16.49% 8.50% 9.71% 172.49% 14394.04% 9.66% e Tab.4.15 indicates that large discrepancies between the two approaches appear mainly in [R]. For [L] the overall differences for the listed combinations of p and r are less e e than 20%. If p , fi , and r are fixed, there is generally a frequency beyond which the differences e e e become larger than a specified tolerance value. This frequency shall be called "threshold frequency." For a 1% tolerance, the threshold frequencies for different p and r are listed e e in Tab.4.16. By converting the threshold frequencies together with the corresponding p , e fi , and r into a parameter r /8 , it is found out that for a given r , this parameter r /6 e e e e e e e Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM is almost constant for different p . With r =24mm, r /S is 0.0213 in [R] and 0.0674 in e e e e [L]. This means that if r /S <0.0213, the differences between Pollaczek's formula and e e FEM are less than 1%. r /8 related to each r is listed at the bottom of Tab.4.16. e e e Table 4.16: Threshold / with Pollaczek's formula for maximum differences <1% Pe 1000 100 10 1 0.1 r./6. r =24mm in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 200kHz > 1MHz 20kHz 200kHz 0.06744 0.02133 s r =100mm in R in L > 1MHz > 1MHz > 1MHz > 1MHz 100kHz > 1MHz 10kHz 200kHz 1kHz 20kHz 0.01987 0.08886 e r,=250m in R in L > 1MHz > 1MHz 100kHz > 1MHz 10kHz 400kHz 1kHz 40kHz 100Hz 4kHz 0.01571 0.09935 r,=1000mm in R in L 60kHz > 1MHz 6kHz 400kHz 600Hz 40kHz 60Hz 4kHz 6Hz 400Hz 0.01539 0.12566 For a 10% tolerance the threshold frequencies are shown in Tab.4.17, and for a 30% tolerance the threshold frequencies are shown in Tab.4.18. Table 4.17: Threshold / with Pollaczek's formula for maximum differences< 10% Pe (fim) 1000 100 10 1 0.1 r /6 t t r,=24rrun in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz — — r =100mm in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 600kHz > 1MHz 60kHz 400kHz 0.15391 0.39738 e r,=250m in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 100kHz > 1MHz 10kHz 100kHz 0.15708 0.49673 r«,=1000mm in R in L > 1MHz > 1MHz 400kHz > 1MHz 40kHz > 1MHz 4kHz > 1MHz 400Hz > 1MHz 0.12566 — Table 4.18: Threshold / with Pollaczek's formula for maximum differences<30% Pe (fim) 1000 100 10 1 0.1 r /6 e e r =24mm in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz e — — r =100mm in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz e — — r =250m in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 800kHz > 1MHz 1.40496 — e r =1000mm in R in L > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 400kHz > 1MHz 40kHz > 1MHz 1.25664 — e From the results presented in this subsection, it can be concluded that Pollaczek's formula for the self impedances of shallowly buried SC cables gives reasonably accurate 92 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 93 results for p in the range of lOOOfim to lfim, and for r in the range of 24mm to 250mm, e e within the frequency range of 1Hz to 1MHz. 4.5.4 [Z] calculations for a cable layout of arbitrary structure It is assumed that the cable in Fig. 3.9(a) is located in an tunnel as shown in Fig. 4.11(a), with /i=1.5m, /i =0.25m, io=0.5m, p —lflm, and p, =p, . The impedances of the cable 2 e e 0 are calculated with the conventional FEM and with the proposed technique. As there is no definite r here, Pollaczek's formula is used with two approximate r of 250mm and e e 500mm. r&=5m is used with the proposed technique and the corresponding FEM mesh at 600kHz is shown in Fig. 4.11(b). (a) cable layout (b) FEM mesh at 600kHz (r =5m) Figure 4.11: The layout of a tunnel installed cable and the FEM mesh at 600kHz fc The impedances found with three approaches are listed in Tab.4.19. "Pollaczekl" represents the results from Pollaczek's formula with r =250mm, and "Pollaczek2" repe resents the results with r =500mm. Tab.4.19 indicates that for small r /8 , Pollaczek's e e e formula gives results which are almost independent of r . The maximum differences bee tween the proposed technique and the conventional FEM are 4.37% in [R] and 0.91% in [L] from 1Hz to 1MHz. The maximum differences between Pollaczek's formula and Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM Table 4.19: [Z] of the tunnel installed cable from three approaches / approach (Hz) 6 60 600 6000 60000 600000 L (mH/km) R (ft/km) conventional proposed Pollaczekl Pollaczek2 conventional proposed Pollaczekl Pollaczek2 conventional proposed Pollaczekl Pollaczek2 conventional proposed Pollaczekl Pollaczek2 conventional proposed Pollaczekl Pollaczek2 conventional proposed Pollaczekl Pollaczek2 #11 #12 #22 0.0447930 0.0448095 0.0448038 0.0448038 0.102843 0.103030 0.103012 0.103005 0.741197 0.745902 0.748330 0.747710 7.39506 7.50986 7.66821 7.61851 61.1415 60.5850 68.8430 65.5297 362.441 355.463 506.977 361.874 0.00597572 0.00599228 0.00599263 0.00599255 0.0611180 0.0613047 0.0613333 0.0613260 0.642955 0.647661 0.649913 0.649293 6.90314 7.01794 7.17553 7.12583 58.3008 57.7442 65.9976 62.6843 353.617 346.639 498.103 353.000 0.420506 0.420523 0.420459 0.420458 0.475638 0.475825 0.475789 0.475781 1.05642 1.06113 1.06332 1.06270 7.22457 7.33937 7.49709 7.44739 58.2408 57.6843 65.9311 62.6177 353.617 346.639 498.103 353.000 ill 2/12 2/22 2.04825 2.05020 2.05154 2.05154 1.81237 1.81420 1.81533 1.81533 1.54492 1.54557 1.54678 1.54680 1.26626 1.26433 1.26433 1.26459 0.97484 0.96679 0.95450 0.95774 0.79508 0.79609 0.72688 0.75557 1.89584 1.89779 1.89906 1.89906 1.66158 1.66341 1.66467 1.66467 1.42010 1.42075 1.42191 1.42193 1.15864 1.15671 1.15670 1.15696 0.88637 0.87832 0.86600 0.86923 0.71164 0.71264 0.64342 0.67211 1.88918 1.89113 1.89241 1.89241 1.65492 1.65675 1.65802 1.65802 1.41345 1.41410 1.41526 1.41528 1.15293 1.15100 1.15099 1.15125 0.88645 0.87840 0.86606 0.86930 0.71164 0.71264 0.64342 0.67211 the conventional FEM are 30.87% in [#] and 14.25% in [L] with r =250mm, and are e 14.88% in [#] and 6.95% in [L] with r =500mm. If p is lOOfhn, the maximum differe e ences become 5.23% in [#] and 0.42% in [L] with r =250mm, and become 4.18% in [#] e and 0.38% in [L] with r =500mm. These results suggest that the impedances of a cable e layout of arbitrary structures can still be calculated with Pollaczek's formula by using an approximate r . e For a loop current of 1-fjO A at 600kHz between the centre conductor of the cable and the earth, the J contours found with the proposed technique and with the conventional FEM are plotted in Fig. 4.12. The plotted area is 2m horizontally by 2.5m vertically. Fig. 4.12 shows good agreement between the field distributions of the proposed technique and of the conventional FEM. It also shows that the tunnel structure does not have a 94 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 95 strong impact on the field distribution, although small deformations can be observed in JR near the tunnel. J .=0.37006-0.01807-(»-l) fl / =0.36736-0.01921-(«-l) iii 32. J =0.83914-0.04282-(«-l) A conventional jr =0.84298-0.04604-(t-l) Ji proposed (r&=5m) Figure 4.12: J distributions in the earth at 600kHz from the FEM 4.6 Summary In this chapter, a technique is proposed to reduce the earth region when the earth penetration depth is large. Good agreements are observed for the calculated impedances and for the field distributions between the conventional FEM and the proposed technique. The solution region is small and independent of frequency, which helps in the meshing Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 96 process. The comparisons show that the proposed technique requires less CPU time than the conventional FEM if partial earth return currents are calculated only once. The conventional FEM and the proposed technique are applied to deeply and shallowly buried cables and tunnel installed cables. In the impedance calculation of deeply buried SC coaxial cables, accurate results can be obtained with the conventional FEM if r > 3»5. With the proposed technique, accurate results are obtained if ri,/8 is small 0 e e at small r , or if r /S is small at large r . For r =24mm, the errors of the proposed e e e e e technique are less than 1% compared with the conventional FEM if r0/t$e < 0.2. In the impedance calculation of shallowly buried and tunnel installed SC coaxial cables, r > 12S is required for the conventional FEM. Comparisons with Pollaczek's b e formula for the shallowly buried cables show that there are discrepancies with the conventional FEM when the earth penetration is small. For typical ranges of p and r , e e however, the discrepancies are reasonably small. With p varying between lOOOJlm to e lfim and with r varying between 24mm to 250mm, the maximum discrepancies are less e than 21% in [R] and less than 9% in [£]. Thefieldsolutions of a tunnel installed cable show that arbitrary structures inside the earth do not have a strong influence on the field distribution. Pollaczek's formula can still be applied tofindimpedances of such a cable by using an approximate r . e Chapter 5 Admittance Calculation with Finite Element Method 5.1 Introduction If the insulating materials in a multiconductor system have complicated geometries, then [Y] of the system cannot be calculated analytically. Instead, numerical methods have to be used. In general, shunt conductances among conductors are ignored. Therefore, the task of the [Y] calculation is simplified into a [C] calculation. According to the assumptions given in Section 2.2, only surface charges exist, and the parallel conductors have uniform cross section longitudinally. This simplifies the capacitance calculation into a two-dimensional steady-state electricfieldsolution problem. The solution region is set up by removing all the regions inside conductors from the solution region in the [Z] calculation. Several numerical methods could be applied, including BEM and FEM. Because the steady-state electric field solution is a special case of the quasi steady-state magnetic field solution, most of the FEM techniques for solving the magneticfieldcan be applied directly to solve the electric field. The corresponding software can easily be adapted to handle the electricfieldsolution. Therefore, only FEM is used in this thesis. In this chapter the general procedures for applying FEM to the steady-state electric field solution arefirstdiscussed. The energy method and the surface charge method to calculate the capacitances from thefieldsolutions are derived next. The general form of [C] for SC coaxial cables is also given. The results of the [C] calculation of a SC coaxial 97 Chapter 5. Admittance Calculation with Finite Element Method 98 cable by FEM show that isoparametric elements are better than simplex elements in this case. 5.2 Principal Equation and F E M solution The assumption that there is no volume charge inside the conductors is justified in [10], based on the equation V . E = £e = i<T v . J = - ia ^dt (5-1) where c is permittivity, p is volume charge density, and a is conductivity. It is pointed out in [10] that, as p = poe ^/^ , any charges introduced into a conductor will dilute - 1 themselves to the surface with a time constant of e/cr. For a poor conducting earth with p = 1000 fim, the corresponding time constant is still a very small value in the order of e 10~ s. In cable related transient studies, the smallest time step is typically in the order 8 of 10 -6 s. Compared to the above time constant, such a time step is large enough to ensure that there is no volume charge inside the conductors. Because the charges are on the conductor surfaces only and there is nofieldinside the conductors, the solution region for the electricfieldwill be composed of insulations bounded partly by conductor surfaces. Introducing the scalar potential <f> as E = -V0 (5.2) the following Laplace equation can be derived eV 0 = 0 2 (5.3) This is the governing equation describing the electricfield.The boundary conditions are still Dirichlet boundary conditions on To and homogeneous Neumann boundary conditions on Y\. Chapter 5. Admittance Calculation with Finite Element Method 99 <f> also appeared in the [Z] calculation. In that case, <f> is the longitudinal potential in conductors which is caused by applied conductor currents. It is constant in the transverse direction in each conductor. In the [C] calculation, however, < > / represents the transverse potential function caused by voltages applied on the conductor surfaces. Comparing the governing equations (5.3) in the [C] calculation with (2.17) in the [Z] calculation, it can be seen that (2.17) will be reduced to the same form as (5.3) if there is no conductor region. With suitable boundary conditions, the real part of the A solution for (2.17) will be the <f> solution for (5.3). That is the reason why the steady-state electric field solution is a special case of the quasi steady-state magnetic field solution. With FEM, <f> is assumed as <f> = £ <f> <p n (5.4) n n=l where 4> is the value of <f> at FEM node n, and NT and <p are the same as before. With n n the same procedure as discussed in Chapter 2, the following algebraic equation can be derived [U}[<f>] = [0] (5.5) where [d>] = U mn = [<j>i,d>2,...,d>NT] T (5.6) f eV<p -V(p ds Js m (m = 1,2,...,N;n = 1,2,...,N ) n T (5.7) R SR and N remain the same as before, and c is the permittivity. By dividing SR into elements and assuming that <f> has the following form in element Ei NBi <f>=Y.€ f ds i Ei in S (5.8) EI n=l where <f> * is the node value of <j> in the element, the integral U E mn U*=e [ Ei Js . B VvZ-Vrf'ds in 5^ becomes (5.9) Chapter 5. Admittance Calculation with Finite Element Method 100 where m — 1,2,..., Ns excluding boundary nodes, n = 1,2,..., JVj^, and i=l,2,... ,M. t and M remain the same as before. Obviously, all. the discussions and formulas in Chapter 2 and Chapter 3 referring to the formulation of [U] for different elements and to the solution of the final algebraic equations are applicable here, provided that i is replaced by c. 5.3 [C] Calculation from the Field Solutions [C\ can easily be calculated from the solved potential field distribution under specific boundary conditions. Similar to the [Z] calculation, there are two methods for finding [C] : the energy method and the surface charge method. Before discussing these two methods, the capacitances of a multiconductor system are first defined: C, is the direct capacitance per unit length between conductor i and the i0 reference conductor, and C mii is the direct mutual capacitance per unit length between conductors i and j. These capacitances are shown in Fig. 5.1. Figure 5.1: Direct capacitances for multiconductor systems Chapter 5. Admittance Calculation with Finite Element Method 101 For conductor i shown in Fig. 5.1, the following equation is derived f ) (juC^dzXVi-Vj^Ii Ii-tiwC.ndzW- + dli (5.10) = l k*i or k= 1 k=1 k?i FC=1 k^i where Cik is the element in [C] , which is related to the direct capacitances by Cu = c, + f ; c i0 (i = i , 2 , . . . , i o (5.12) (i,k = l,2,...,K;k^i) (5.13) mii fc^i C ik = -C mih For the complete system, (5.11) has the matrix form This equation is essentially (2.2) provided that [G] is ignored. (5.12) and (5.13) are used in the energy method to find [C] from thefieldsolutions. 5.3.1 The energy method With this method, the electric energy stored in thefieldunder specific boundary conditions is calculated, and the elements of [C] are derived from the energy. Once the potential distribution is known, the following equation is used tofindthe electric energy from the field = \ I cE • Eds = i / i = eV<p • V<pds M (5-15) Chapter 5. Admittance Calculation with Finite Element Method 102 where [^] = (5-16) [<j> ] is the node value vector of <f> in element Ei. c^. is the permittivity in Ei. Ei Fig. 5.2 shows the cross section of the conductor system in Fig. 5.1. From circuit analysis, if Vi = Vo ^ 0 and Vj = 0 (j = 1,2,..., K; j ^ i), the stored electric energy is Figure 5.2: Direct capacitances under DC condition given by (5.17) k= k*i i For these conductor voltage conditions, the corresponding boundary conditions in the field solution are that all potentials at the boundary nodes on the surface of conductor i are Vo and that the potentials at the rest of the boundary nodes are zero. If the energy calculated from the potential distribution is represented by WJf , Cn becomes c By assuming Vi = Vj = V ^ 0 and 0 2W ^ (5.18) = 0 (k = 1,2,..., K; k ^ i; k ^ j), the stored Chapter 5. Admittance Calculation with Finite Element Method 103 electric energy from circuit analysis becomes wl = \c ,X+ - £ c ,X-\c v 1 s m mtJ 2 0 k=i 1 1 2 1 K 2 2 fc = 1 = \(C ii + C + 2C )V jj ij (5.19) 2 0 For the field solution, the potentials at the boundary nodes on the surfaces of conductors i and are Vo- The potentials at the rest of the boundary nodes are zero. If represents the stored energy in the field, C*j will be ~ 14? 2 Y ( 5 - 2 0 ) As the system is linear, the field needs to be solved only K times, with one conductor surface having non-zero potential each time. From these K solutions, Wg is calculated F and Cn is evaluated. By superimposing two solutions, W Ep and C can then be found. t J This method is simple, and the potential distribution solved by F E M is used directly to find [C] . One interpretation of F E M for the Laplace equation is that F E M trys to find a solution that will minimize the energy in the field. Therefore, [C] calculated with the energy method should always be slightly larger than the exact value, unless the solution itself is exact. 5.3.2 The surface charge method With this method, [C] is found from the surface charges per unit length on the conductors under specific boundary conditions. The surface charges are calculated from the integration of D along the contours of the conductors. Chapter 5. Admittance Calculation with Finite Element Method 104 From electrostatic field analysis, the surface charges per unit length on the conductors in Fig. 5.2 are related to [C] through the following equation [q] = [C][V] (5.21) where [<z] = [<Zi,<Z,...,<Zx] (5-22) r 2 qi is the surface charge per unit length on conductor i. [V] is the same as before. If Vi = VQ ^ 0 and V - = 0 (j = 1,2,..., K;j ^ i), the elements of the t'th column in [C] } can be found from the corresponding surface charges as C* = # (j = l,2,...,iiQ (5.23) qj is given by q = f } D • dl = / ±\D\y/dx + 2 dy 2 (5.24) where YCJ is the periphery of the cross-section area of the jth conductor, dl is an integral element with a direction normal to Tc The plus and minus signs in the above formula r are related to the direction of D. ID I is calculated from (»Y (»)'+ | D | - 4 E | - « | V * | - « ^ = j+ 3 ^ (5.25) The system of equations needs to be solved K times, with one conductor surface having non-zero potential each time. In order to evaluate the integral (5.24) in a simplex element, the vertices are numbered in the same way as shown in Fig. 3.1. Three possible combinations can be obtained by permuting three node numbers: (1,2,3), (2,3,1), and (3,1,2). (m,n, o) is used to represent these three combinations. Chapter 5. Admittance Calculation with Finite Element Method 105 Assuming that the integral (5.24) in element Ei is along the side from vertices m to n on which £<, = 0 and d(, = 0, the corresponding integral Aq can be derived as 0 A<? = e yJ(x -x )* Ei n + (y -y y n m ±^ (^j £ m + (^j (5.26) d( n and dx ~du dx - — d u — - P (N ,( ) jm p £. > \dU 9 P i n { N dx M ( 5 — m — where P , P , and P are given by (3.5) and (3.6). jm jn jo also be expressed in terms of and SC„ dx. + - 2 9 ) (5.30) and in (5-27) and (5.28) can or in terms of and if (5.29) and (5.30) are modified correspondingly. In order to simplify the calculation of -g^- and P (N ,Q m p is changed into real polynomial form Pm(N ,0 p 1 m (5.31) = —52a C "•m i mi =0 in which d is the common denominator. a - is the polynomial coefficient. d and a m mt m mt are associated with order N and with m, and they are listed for order 1 to 6 in Tab.5.1. p As | £ and | ^ can be found for any given (( , £„, ( ) from (5.27) to (5.30), Aq in (5.26) m 0 can now be evaluated by applying the Gaussian quadrature discussed in Section 3.3.1. Chapter 5. Admittance Calculation with Finite Element Method Table 5.1: d and a m 1 2 3 4 5 6 in P (N ,Q mi m dm 0 l 1 1 1 0 0 l 2 0 1 2 3 0 1 2 3 4 0 1 2 3 4 5 0 1 2 3 4 5 6 1 1 1 1 1 2 2 1 1 1 3 3 1 1 2 6 24 24 1 1 1 1 2 5 10 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 m Oml 3 -3 2 9 -9 9 4 -2 4 -3 8 -24 22 32 -48 32 5 -5 10 -30 24 25 -75 275 -250 125 -750 875 625 -1250 625 6 -3 2 -3 6 -10 18 -18 33 -75 137 36 -108 315 -675 108 -540 1530 324 -1620 dy Sj/ d£ du §± L dy J Om6 d_ dx 648 and | £ in element E{ can be derived as 8y dv r <&m5 2 dx dx dv dv —- «m4 2 -1 dx dx Ora3 1 For isoparametric elements, dx, dy, §± =.£Efeo«-m<C* p <*m2 106 - dv du (5.32) (5.33) L where [J ] and [D ] are given by (3.22) and (3.25), respectively. \0\ is given by (3.19) for Ei Ei quadrilateral isoparametric elements, or by (3.28) for triangular isoparametric elements. If the integration within a quadrilateral element is along the side from local nodes 1 Chapter 5. Admittance Calculation with Finite Element Method 107 to 2 as shown in Fig. 3.5, u = — 1 and du = 0, and from (5.32) dx dv 1 in 0 dv (5.34) The corresponding integral Aq becomes (5.35) Once the local coordinates of sampling points for the numerical integration are known along this side, | £ , jjj, and jjjj are calculated from (5.33) and (5.34), respectively, and Aq from (5.35). For integrations along the other sides, similar formulas can be derived. If the integration is within a triangular element along the side from local nodes 1 to 2 as shown in Fig. 3.6, u = — v + 1; therefore, dx dv ±L J* dv J with A 5 having the same form as in i du L dv (5.36) - i (5.35). The concept of this method is very simple; however, the software implementation is rather complicated compared with the energy method. Although <f> is continuous in FEM solutions, E is not. This may introduce errors in [C] . 5.4 [C] Calculation of SC Coaxial Cables In this section the capacitances of SC coaxial cables are obtained, and the results are compared with those from analytical formulas. 5.4.1 General form of [C] for SC coaxial cables The same notations will be used as in Section 2.6 or as shown in Fig. 2.4. For the SC cable shown in Fig. 2.4 there are K insulations. Thefcthinsulation is between conductors Chapter 5. Admittance Calculation with Finite Element Method 108 k and k+l, and the capacitance per unit length related to the insulation is 2«fc (A? = 1,2,..., K) (5.37) where e is the permittivity of thefcthinsulation. The general form of [C] for SC coaxial k cables, which only has three diagonals[10], is c{ cf cf ci riod °3 [C} = 0 Cf rid °3 riod °4 (5.38) riod rid riod riod rid where cf = C>IN, C/Jv _! + CiN (k = 2,3,..., K) —CiNu.x (k — 2,3,..., if) t c? 5.4.2 k (5.39) [C] calculation of a SC coaxial cable by F E M [C] of the SC coaxial cable shown in Fig. 3.9 is calculated by FEM. As shown in the figure, e,=l is assumed for both insulations. The solution region is similar to Fig. 3.9(b), except that the two conductor regions are removed. As in the [Z] calculation, the span angle 0 and the division radii will affect the results. When studying the influence of division radii on [C] , a small 0 is used. In Tab.5.2 five different radius divisions are listed. As the solution region is made up by two disconnected insulation regions, the division radii for each division scheme are listed for these two regions separately. Chapter 5. Admittance Calculation with Finite Element Method 109 Table 5.2: Radius divisions in [C] calculation of the SC coaxial cable division division 1 division2 division3 division4 division5 12 12 12 12 12 division radii (mm) inner insulation outer insulation 18 22 24 22 23 24 15 18 14 16 18 22 22.7 23.3 24 13 15 17 18 22 22.3 23 23.7 22 22.1 23.9 24 12.3 17.7 18 24 The results of the [C] calculation are given in Tab.5.3 for different division schemes with 6 = 1°. Values from the analytical formula (5.38) are also included in the table. It can be seen that, with all division schemes, the capacitances calculated by the energy method have four digit accuracy for all types of elements except the first order simplex element. With this method, the calculated values for C\i and C i are the same as — C . 2 u The capacitances calculated by the surface charge method depend on division scheme. C\2 is the same as — C n and C i is slightly different from C12 due to different integration 2 surfaces in the surface charge calculation. It can be seen from Tab.5.3 that for simplex elements higher order methods give more accurate results, at the expense of more nodes and of longer computation time. For all the division schemes, the 3rd order simplex element seems to be an optimum choice, since the results are not improved significantly by further increasing the order of the element. For the numerical integration in the surface charge method the number of sampling points in the Gaussian quadrature can easily be decided for a simplex element. For a N th order simplex element, with <j> being a polynomial function of order N , E calculated p p from <f> by differentiation will be a polynomial function of order N — 1. Therefore, the p number of sampling points will be ((N + l)/2), truncated to the nearest integer. For p isoparametric elements, the results do not change very much after the number of sampling points goes beyond five. Tab.5.3 shows that different radius division schemes will mostly affect the results Chapter 5. Admittance Calculation with Finite Element Method 110 Table 5.3: [C] of a two-conductor SC coaxial cable element analytical ISO divisionl sim2 sim3 sim4 sim5 sim6 ISO division2 siml sim2 sim3 sim4 sim5 sim6 ISO division3 siml sim2 sim3 sim4 sim5 sim6 ISO division4 siml sim2 sim3 sim4 sim5 sim6 division5 siml sim2 sim3 sim4 sim5 sim6 ISO Energy Method (/iF/km) Surface charge method (/iF/km) Cu (-C121 -C21) Cn (-C12) Cn 0.137017 0.775503 0.137017 -0.137017 0.775503 0.137037 0.137039 0.137018 0.137017 0.137017 0.137017 0.137018 0.137503 0.137019 0.137017 0.137017 0.137017 0.137017 0.137017 0.137236 0.137018 0.137017 0.137017 0.137017 0.137017 0.137017 0.137178 0.137017 0.137017 0.137017 0.137017 0.137017 0.137029 0.138374 0.137029 0.137017 0.137017 0.137017 0.137017 0.775523 0.775533 0.775507 0.775505 0.775504 0.775504 0.775504 0.776107 0.775509 0.775505 0.775504 0.775504 0.775504 0.775503 0.775784 0.775507 0.775505 0.775504 0.775504 0.775504 0.775503 0.775717 0.775505 0.775504 0.775504 0.775504 0.775504 0.775515 0.777170 0.775517 0.775504 0.775504 0.775504 0.775504 0.133333 0.133337 0.136559 0.136963 0.136999 0.137014 0.135886 0.122225 0.135882 0.136932 0.137001 0.137004 0.137004 0.136476 0.126680 0.136467 0.136978 0.137000 0.137000 0.137000 0.136871 0.131691 0.136852 0.136987 0.136988 0.136987 0.136989 0.137015 0.136666 0.136955 0.136938 0.136948 0.136966 0.136981 -0.133333 -0.133337 -0.137476 -0.136963 -0.137035 -0.137014 -0.136261 -0.150004 -0.136274 -0.137069 -0.137026 -0.137029 -0.137029 -0.136701 -0.145309 -0.136719 -0.137045 -0.137035 -0.137036 -0.137036 -0.136942 -0.141097 -0.136976 -0.137053 -0.137053 -0.137052 -0.137048 -0.137022 -0.139537 -0.137140 -0.137119 -0.137091 -0.137064 -0.137049 0.771014 0.771034 0.775724 0.775439 0.775278 0.775487 0.774537 0.774416 0.774382 0.775334 0.775290 0.775312 0.775346 0.775082 0.773857 0.774854 0.775232 0.775239 0.775291 0.775349 0.775409 0.775312 0.774867 0.774995 0.775182 0.775321 0.775380 0.775506 0.776884 0.774046 0.775028 0.775356 0.775324 0.775392 Chapter 5. Admittance Calculation with Finite Element Method 111 from the surface charge method. If the potentials of the centre conductor and on the outmost boundary in Fig. 3.9(b) are zero, and the potential of the second conductor is 1 V, the corresponding surface |E| calculated by FEM with isoparametric elements are given in Tab.5.4. The table shows that a closer division radius towards the boundary Table 5.4: Surface |E| for different divisions with isoparametric elements analytical divisionl division2 division3 division4 division5 r=12 mm 0.205525 0.200000 0.203829 0.204713 0.205306 0.205522 |E(r)| (V/mm) r=18 mm r=22 mm r=24 mm 0.137017 0.522398 0.478865 0.133333 0.478261 0.521739 0.478720 0.136261 0.522226 0.522312 0.136701 0.478795 0.478852 0.136943 0.522381 0.478863 0.137023 0.522396 will improve the surface E results and consequently improve the accuracy of the surface charge method, as indicated in Tab .5.3. When 9 varies, the FEM mesh will be changed. If the whole circular region in Fig. 3.9(a) is considered, different 9 will create different numbers of nodes and numbers of elements. The relative errors of numerical results compared with those given by the analytical formula are plotted in Fig. 5.3 for isoparametric elements and for the 3rd order simplex element. In thefigure,EM stands for the energy method and SCM for the surface charge method, iso and sim3 are the same as before. It can be seen that isoparametric elements produce accurate results even at very large span angles, while the 3rd order simplex element only produces accurate results at small span angles. For other orders the results are similar. Therefore, isoparametric elements are more computation efficient, as less computation time will be required by fewer nodes and fewer elements in the FEM mesh. Chapter 5. Admittance Calculation with Finite Element Method 0.8 0.4 3! v = C n (divisions) s 20 —r 60 400 (degree) 120 50 5 40 d 60 EM sim3 - n = C n 7division3) •v = C u (division5) = C M (division3) -x = C u (division5) 50£i 40 d 30- 3 0 CM CM 60 BO 120 (degree) SCM sim3 -a = C u (diviaion3j v = C u (division5) ••o = C»i (division3) -+ = Cti (division5j •A = C M (division3) -x = C M (division5) 20 u o H 0.0 = = 70 80 H • X 9 70 1 H= A d 0.4-fe & A 0.2 0.1 0.0 o.e- s -A = C M idiviaion3) x = C M (division5) d 0.2 oU n SCM_iso division3 — Oil (division5) -» = Cn division3, division5 Cn division3. -D = C u EM_iso a = C u (division3) 112 10 H 20 w 20 40 10 8 eo eo (degree) 100 120 0 8 (degree) 120 0 5.3: Errors in [C] calculation of a SC coaxial cable for different span angles Figure 5.5 Summary In this chapter general procedures for solving the electrostatic field with FEM are discussed. The energy method and the surface charge method for calculating [C] from the field solutions are derived. The energy method is simple and easy to implement compared with the surface charge method. Both methods are applied to the SC cable shown in Fig. 3.9. For this case the results show that isoparametric elements achieve higher accuracy with fewer elements and nodes in the FEM mesh than simplex elements. The results also show that the Chapter 5. Admittance Calculation with Finite Element Method 113 energy method has less stringent requirements on the mesh, and has higher accuracy than the surface charge method. For isoparametric elements a division radius close to the conductor surface will improve the [C] results found by the surface charge method, while for simplex elements similar improvements can be achieved by increasing the radial divisions. Chapter 6 Case Studies in [Z] and [C] Calculations 6.1 Introduction There are analytical formulas for the calculation of parameters of most types of power cables. These formulas are generally derived with certain approximations. For example, the impedance formulas of PT cables are derived by replacing the conductors inside the pipe with current filaments located at the centres of the conductors. By applying FEM, the parameters can be calculated without some of these approximations. This provides a way of verifying the validity of analytical formulas. In this chapter, FEM is applied to the [Z] calculation of buried or tunnel installed multiphase SC cables, of PT cables, of sector-shaped cables, and of stranded conductors. Capacitances of PT cables and sector-shaped cables are also calculated with FEM. The numerical results of [Z] and [C] are compared with those from the analytical formulas. The comparisons show that for buried multiphase SC cables, accurate self and mutual impedances can be obtained with Pollaczek's formula. For tunnel installed SC cables, reasonably accurate self and mutual impedances can still be obtained with Pollaczek's formula. With PT cables, the analytical formulas may introduce errors at high frequencies due to neglected proximity effects. For sector-shaped cables, FEM results are compared with those from an approximate analytical formula recently suggested by Ametani[40]. Discrepancies are observed. For stranded conductors, comparisons are made among the results found with FEM, the subdivision method, the GSW formula, and Borges da Silva's 114 Chapter 6. Case Studies in [Z] and [C] Calculations 115 formula. These comparisons show that close agreement is obtained between FEM results and the GSW formula if a factor of 1.6 is used in the formula instead of 2.25. Good agreement exists between FEM results and Borges da Silva's formula. 6.2 [Z] Calculations of Buried or Tunnel Installed Multiphase Cable Systems In Section 4.5 the impedances of a single phase cable are calculated, where mutual impedances only exist between the conductors within the cable. For multiphase cable systems, the mutual impedances among the conductors of different phases must be calculated as well. To use the formula (4.6) for the impedance between two phases, (0,h) and (xp,yp) are assumed to correspond to the locations of the two phases, respectively. In this section, the impedances of two 230kV three-phase SC cable systems are calculated with (4.6), with the conventional FEM, and with the proposed technique. Each phase consists of a two-conductor SC cable, as shown in Fig. 6.1(a). One system is buried as shown in Fig. 6.1(b), and the other is installed in a tunnel as shown in Fig. 6.1(c). The FEM meshes around the cables at 60Hz are given in Fig. 6.2. As there are six conductors sheath: 0 = 4800 S/mm core: o=57000 S/mm •s r 1 2 (a) cable data earth earth u ^i=12.1mm '« =24.5mm fr2=42.3mm B=46mm rA3=50.9mm air air r, \ s r e A B C © © 1 h=l.2m j=300mm r=50.9mm p.= lOOQm LU=u ^ e e © 0 n A=1.2m h =400mm /j=100mm S =300mm w =150mm 2 3 p =100£2m \l =\i g e Q (c) tunnel installed system (b) buried system Figure 6.1: 230kV three-phase cable systems Chapter 6. Case Studies in [Z] and [C] Calculations 116 (a) buried cable system (b) tunnel installed cable system Figure 6.2: Meshes around the cables at 60Hz for the two systems in each system, the final [Z] has the form of [ AA] [ AB1 [ZAB] [ * B B ] [*BC] Z 0]= . i AC\ Z [ AC1 Z i BC\ Z Z t-) 6 1 [ CC\ Z Subscripts A, B, and C represent the phases. All submatrices are 2x2 matrices, and [•^AA1» [^BB1» an< ^ [^CCl a r e syn^etric. The sheath conductor is numbered after the core conductor in each phase. The impedances of the buried cable system are listed in Tab.6.1 and Tab.6.2. Because the meshes are symmetrical with respect to the center cable, [^AA] ^CCl = anc * [Z^g]=[ZgQ]. With Pollaczek's formula, it is assumed that [^AA] [^BB] [^CC3> * * = all four elements in [^AB! a r e same > ana ^ that all f ° u r elements in [^ACl = a r e na ^ n e s a m e - Tab.6.1 and Tab.6.2 confirm that these assumptions are reasonably accurate for the given p =100flm. The tables show that the differences between the results from Pollaczek's e formula and from FEM are small at low frequencies, but noticeable at high frequencies. Good agreement is observed between the conventional FEM and the proposed technique. The maximum differences between Pollaczek's formula and the Chapter 6. Case Studies in [Z] and [C] Calculations Table 6.1: [•ZAAKI^CC'D [^BBI 117 °^ *^e buried three-phase cable system L (mH/km) R (ft/km) / Rn (Hz) 60 6000 600000 [*AA1 conv. prop. conv. [*BB! prop. Pollaczek [*AA! conv. prop. [*BB1 conv. prop. Pollaczek [*AA! conv. prop. t*BBl conv. prop. Pollaczek R12 0.0748300 0.0749538 0.0750753 0.0751972 0.0745183 6.37885 6.39170 6.38720 6.39999 6.37657 694.949 696.979 694.819 696.919 698.287 Table 6.2: [ ^ A B K I ^ B C J ) R22 0.0597006 0.262705 0.0598244 0.262829 0.0599458 0.262960 0.0600678 0.263072 0.0594093 0.262382 6.16838 6.33678 6.18123 6.34963 6.17672 6.34512 6.18952 6.35792 6.16547 6.33386 691.010 691.011 693.041 693.041 690.880 690.881 692.980 692.981 694.327 694.327 t^ACl °^ £11 L12 2.10829 2.11098 2.10904 2.11172 2.12015 1.61213 1.61402 1.61034 1.61225 1.62790 1.09676 1.09769 1.09482 1.09576 1.11182 1.96073 1.96342 1.96147 1.96416 1.97267 1.49146 1.49335 1.48967 1.49158 1.50717 0.98654 0.98747 0.98460 0.98554 1.00154 buried three-phase cable system L (mH/km) J2 (ft/km) / (Hz) [*AB1 60 f*ACl [*AB1 6000 f*ACl PABI 600000 [*ACl conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. Ru 0.05944 0.05957 0.05940 0.05895 0.05908 0.05940 6.085 6.098 6.092 6.066 6.079 6.091 685.9 687.9 689.1 679.9 681.8 682.5 R12 -R21 L22 1.95793 1.96062 1.95868 1.96137 1.96988 1.48897 1.49086 1.48718 1.48909 1.50468 0.98654 0.98747 0.98460 0.98554 1.00154 R22 0.05944 0.05957 0.05944 0.05957 0.05944 0.05957 0.05895 0.05908 0.05895 0.05908 0.05895 6.085 6.098 6.085 6.098 6.085 6.098 0.05908 6.066 6.066 6.066 6.079 6.079 6.079 685.9 687.9 685.9 687.9 685.9 687.9 679.9 681.8 679.9 681.8 679.9 681.8 ill Ll2 L21 L22 1.586 1.589 1.589 1.449 1.452 1.451 1.120 1.122 1.125 0.989 0.990 0.986 0.622 0.623 0.626 0.492 0.492 0.488 1.586 1.589 1.586 1.586 1.589 1.589 1.449 1.452 1.449 1.449 1.452 1.452 1.120 1.122 1.120 1.120 1.122 1.122 0.989 0.990 0.989 0.990 0.989 0.622 0.623 0.622 0.623 0.622 0.623 0.492 0.492 0.492 0.492 0.492 0.492 0.990 conventional FEM are 1.8% in [R] and 1.78% in [L] from 1Hz to 1MHz, and the maximum differences between the proposed technique and the conventional FEM are 0.49% in [R] and 0.18% in [L] from 1Hz to 1MHz. The impedances of the tunnel installed cable system are listed in Tab.6.3 and Chapter 6. Tab.6.4. Case Studies in [Z] and [C] Calculations 118 r =300mm is used with Pollaczek's formula. Again, the differences between e Pollaczek's formula and the conventional FEM become noticeable at high frequencies. Table 6.3: [^AAK^CC]) a n < * I^BBI u n n e l installed three-phase cable system L (mH/km) # (ft/km) / (Hz) Rll 60 6000 600000 Table 6.4: °^ ^e * 1*AA1 conv. prop. [*BB1 conv. prop. Pollaczek 1*AA1 conv. prop. [*BB1 conv. prop. Pollaczek [*AA1 conv. prop. conv. [*BB! prop. Pollaczek [-^ABKI^BC^) 0.0748252 0.0749517 0.0750711 0.0751999 0.0745182 6.37058 6.38852 6.37804 6.39624 6.37618 675.093 685.630 672.579 683.822 695.484 a n c #12 #22 0.0596958 0.0598223 0.0599417 0.0600705 0.0594092 6.16010 6.17804 6.16757 6.18576 6.16508 671.155 681.692 668.641 679.884 691.523 0.262700 0.262827 0.262946 0.263075 0.262382 6.32850 6.34644 6.33597 6.35416 6.33347 671.155 681.692 668.641 679.884 691.523 ^ l^ACl °^ L\\ 2.11131 2.11394 2.11116 2.11376 2.12015 1.61540 1.61716 1.61276 1.61448 1.62790 1.10341 1.10176 1.10084 1.09903 1.11193 PABI 60 [*ACl I*AB1 6000 [*ACl [*AB! 600000 t^ACl conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. conv. prop. Poll. L22 tunnel installed three-phase cable system L (mH/km) # (ft/km) / (Hz) Lu 1.96375 1.96096 1.96637 1.96358 1.96360 1.96080 1.96620 1.96340 1.97267 1.96988 1.49473 1.49224 1.49649 1.49400 1.49209 1.48360 1.49381 1.49132 1.50717 1.50468 0.99319 0.99319 0.99154 0.99154 0.99061 0.99061 0.98881 0.98881 1.00166 1.00166 #11 #12 0.05944 0.05957 0.05940 0.05894 0.05907 0.05940 6.077 6.095 6.092 6.060 6.078 6.091 665.5 676.0 689.1 661.4 671.9 682.5 0.05944 0.05957 0.05944 0.05944 0.05957 0.05957 #21 #22 0.05894 0.05907 0.05894 0.05804 0.05907 0.05907 6.077 6.095 6.077 6.095 6.077 6.095 6.060 6.078 6.060 6.078 6.060 6.078 665.5 676.0 665.5 676.0 665.5 676.0 661.4 671.9 661.4 671.9 661.4 671.9 Lu L12 L21 L22 1.584 1.587 1.589 1.584 1.584 1.584 1.587 1.587 1.587 1.450 1.451 1.118 1.120 1.125 0.987 0.989 0.986 0.624 0.622 0.626 0.494 0.492 0.488 1.450 1.450 1.450 1.448 1.448 1.448 1.448 1.118 1.118 1.118 1.120 1.120 1.120 0.987 0.989 0.087 0.989 0.987 0.989 0.624 0.622 0.624 0.622 0.624 0.622 0.494 0.492 0.494 0.492 0.494 0.492 Chapter 6. Case Studies in [Z] and [C] Calculations 119 The maximum differences are 4.81% in [R] and 1.87% in [L] from 1Hz to 1MHz. The maximum differences between the proposed technique and the conventional FEM are 2% in [R] and 0.75% in [L]. For both systems, the J distributions in the earth from the proposed technique and the conventional FEM are given in Fig. 6.3. In this distribution, a loop current of l+j'O A is assumed between the left phase and the earth at 600kHz. The field plotting area is 3m vertically by 3m horizontally. Comparing the fields of the buried cable system with the tunnel installed cable system shows that the tunnel structure does not have a significant influence on the field. Fig. 6.3 also indicates a good agreement between the fields from the proposed technique and the fields from the conventional FEM. 6.3 [Z] and [C] Calculations of P T Cables 6.3.1 The [Z] calculation of P T cables For the three-phase PT cable with SC coaxial cables inside, as shown in Fig. 1.1(a), it is very difficult, if not impossible, to obtain an analytical field solution. H the conductors inside the pipe are replaced by current filaments, however, the fields inside the pipe and within the pipe conductor can be solved analytically, and approximate impedance formulas can then be derived from the field solutions. The fields inside the pipe and within the pipe conductor were solved first by Tegopoulos et alin 1971[8], and the approximate formulas for the [Z] calculation were developed first by Brown et al in 1976[13]. These formulas are widely used in cable parameter programs [23]. Several of the assumptions in the approximate formulas ignore factors which may influence the results of the [Z] calculation. It is assumed that the fields within each SC coaxial cable inside the pipe remain cylindrical, to be able to use the formulas discussed in Section 2.6. This ignores the influence of the pipe and adjacent cables on the field J = 6903.9 - 47.4-0'-1) Ri , =36453-1340«0-1) hi = 6928.1 - 47.2 < i -1) =36486 - 1339'(j-1) /«, =6760.5 - 44.3 •(/-!) //, =31626-1105-(. -1) A* = 6857.8 - 41.7.(1- / / , =31563-1104'( i-l) conventional proposed conventional proposed (a) buried cable system (b) tunnel installed cable system Figure 6.37 /distributions in the earth at 600 kHz with the earth current being 1+; 0 A Chapter 6. Case Studies in [Z] and [C] Calculations 121 distribution of the concerned cable. As thefieldsolutions were derived for a loop current between a single currentfilamentand the pipe, the influence of the other non-current carrying conductors on thefielddistributions inside the pipe and within the pipe conductor is ignored. The influence of thefinitedimension of the current carrying conductor, which is replaced by a currentfilamentfor the formula derivation, is also ignored. With FEM, these factors can be considered, and their effects of can be studied. A three-phase 230kV PT cable is given in Fig. 6.4. Each phase consists of a two-conductor SC cable as shown in Fig. 6.4(a). The non-linearity of the steel pipe is ignored, and a 240 mm (a) cable data (b) triangle arrangement (c) cradle arrangement Figure 6.4: A 230kV PT cable system constant /*,=500 is assumed instead. Two arrangements of the SC cables inside the pipe are analyzed: the triangle arrangement and the cradle arrangement shown in Fig. 6.4(b) and (c), respectively. The meshes at 6kHz for both arrangements are plotted in Fig. 6.5. In the calculation, the pipe will be used as the return path, and the earth is not included. A=0 is assumed on the boundary in the meshes in Fig. 6.5. The structure of [Z] is the same as in (6.1). For the triangle arrangement the impedances are listed in Tab.6.5 and Tab.6.6. Due to symmetry, [^BB] f^CC^ = [^AB] ^ACl* " = a n a " m *he tables stands for the results Chapter 6. Case Studies in [Z] and [C] Calculations 122 (a) triangle arrangement (b) cradle arrangement Figure 6.5: Meshes for the PT cables at 6kHz found with the analytical formulas [23]. Large differences are observed between the results from FEM and those formulas at high frequencies. The maximum differences are 88% in [R] and 68% in [L] from 1Hz to 1MHz. As mentioned earlier in this section, the formulas ignore the influence of non-current carrying conductors on the field distributions inside the pipe and within the pipe conductor. To study this influence, a 1+j'O A loop current is assumed between the sheath of the upper SC cable and the pipe. \A\ distributions caused by the loop current with or without two lower cables are plotted for 60Hz and 6kHz in Fig. 6.6. The corresponding | J\ distributions on the inner surface of the pipe are given in Fig. 6.7. The A contour lines represent the magneticfluxlines[34]. It can be seen in Fig. 6.6 that the distortion of the |.4.| distribution caused by the presence of the two lower cables is severe at 6kHz but not severe at 60Hz. Similarly, the | J | distribution at 6kHz on the inner surface of the pipe is altered greatly near the two lower cables if the two cables are present, while at 60Hz the \ J\ distribution is only slightly affected by the lower cables. As reflected in the impedances in Tab.6.5 and Tab.6.6, the differences between the results from FEM and those from the analytical formulas are small at 60Hz (< 3.1%) and large at 6kHz (< 26.9%). Chapter 6. Case Studies in [Z] and [C] Calculations Table 6.5: [ - ^ A I A N < ^ [^BBK^CC"]) (Hz) #11 I*AA1 60 [%B! 1*AA1 600 [ZBBI PAAI 6000 13BB] [*AA! 60000 ( BB1 Z I*AAJ 600000 I^BB) FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. L (mH/km) #12 #22 I'll L\2 2/22 0.86707 0.89132 0.90023 0.92884 0.86545 0.88969 0.89862 0.92721 0.260819 0.254153 0.264551 0.257818 0.662920 0.656159 0.666666 0.659823 0.975600 0.908456 1.01266 0.977012 0.916694 0.849402 0.053773 0.917958 1.31873 1.25125 1.35584 1.31981 0.48853 0.53731 0.50442 0.55370 0.36789 0.41660 0.38378 0.43298 0.36628 0.41498 0.38216 0.43137 3.20001 3.20516 3.53022 3.89726 3.00974 3.01413 3.33996 3.70623 3.40539 3.40962 3.73563 4.10172 0.29346 0.36005 0.29957 0.34546 0.18268 0.24922 0.18879 0.23464 0.18108 0.24762 0.18719 0.23304 11.9585 12.4560 16.3469 16.3205 10.6357 11.1361 15.0241 15.0007 10.7013 11.2030 15.0897 15.0675 0.23325 0.29353 0.22250 0.25717 0.12740 0.18764 0.11665 0.15128 0.12666 0.18691 0.11591 0.15055 42.6742 41.3592 76.1295 56.6433 38.1846 36.8655 71.6390 52.1496 38.1927 36.8741 71.6480 52.1582 0.20765 0.26805 0.18052 0.22257 0.10516 0.16552 0.07803 0.12004 0.10516 0.16552 0.07803 0.12004 [^BCl °^ ^ PT cable with a triangle arrangement e L (mH/km) # (ft/km) #11 #12 #21 #22 Lu 2/12 2/2i 2) 2 FEM ana. 0.2498 0.2493 0.2498 0.2493 0.2470 0.2482 0.6858 0.6824 0.6396 0.6857 0.6824 0.6470 0.6396 0.6857 0.6824 0.2482 0.2498 0.2493 0.2470 0.2482 0.6857 0.6824 ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. 0.2498 0.2493 0.2470 0.2482 0.6396 0.6396 0.7763 0.7611 0.7529 0.7425 0.7764 0.7611 0.7528 0.7425 0.7763 0.7611 0.7528 0.7425 0.7764 0.7611 0.7528 0.7425 0.2357 0.2357 0.2055 0.1988 0.2357 0.2357 0.2055 0.1988 0.2357 0.2357 0.2055 0.1988 0.2357 0.2357 0.2055 0.1988 2.506 2.184 2.372 2.010 2.506 2.184 2.372 2.010 2.506 2.506 2.184 2.184 2.372 2.372 2 . 0 1 0 2.010 0.0947 0.1069 0.0722 0.0782 0.0947 0.1069 0.0722 0.0782 0.0947 0.1069 0.0722 0.0782 0.0947 0.1069 0.0722 0.0782 8.729 6.283 7.696 5.428 8.729 6.283 7.696 5.428 8.729 6.283 7.696 5.428 8.729 6.283 7.696 5.428 0.0466 0.0711 0.0282 0.0474 0.0466 0.0711 0.0282 0.0474 0.0466 0.0711 0.0282 0.0474 0.0466 0.0711 0.0282 0.0474 33.58 18.98 33.58 18.98 33.58 18.98 0.0283 0.0604 0.0283 0.0604 0.0283 0.0604 0.0283 0.0604 25.02 25.02 25.02 33.58 18.98 15.86 15.86 15.86 15.86 0.0386 0.0386 0.0386 0.0386 (Hz) I*AB1 [*BCl I*AB1 [*BC] 1*AB1 [*BCl [*AB! 60000 e 0.282226 0.275526 0.285957 0.279190 / 6000 ^ with a triangle arrangement 1.01507 1.03926 1.04822 1.07678 Table 6.6: [ ^ A B K I ^ A C " ] ) 600 c a D # (n/km) / 60 °^ [*BCl 1*AB1 600000 [*BCl 123 FEM 0.2470 25.02 0.6470 2 0.6470 0.6470 0.0136 0.0136 0.0136 0.0136 Chapter 6. Case Studies in [Z] and [C] Calculations 124 one cable three cables |4| = 37.832 • (*-l) = 10.066 • («-l) (a) |A| at 60Hz (/zsV/km) (b) at 6kHz (/xsV/km) Figure 6.6: \A\ distributions caused by the loop current at 60Hz and 6kHz -I 17.0 0.030 -° = one cable v = three cables 16.0 H I -o = one cable 0.025 15.0 & 0.020 14.04 o.oisH 13.0 0.010H 12.0 —i i 1 1 1 1 1— 4S 80 135 180 225 0 (degrees) 270 (a) surface \ J\ at 60Hz 315 380 0.005 0 —i 45 1 <t> 1 1(degrees) 1 1 1— 80 135 180 225 270 315 (b) surface \ J\ at 6kHz Figure 6.7: \J\ on the inner surface of the pipe caused by the loop current 380 Chapter 6. Case Studies in [Z] and [C] Calculations 125 If only the upper cable is in the pipe, the maximum differences between the results from FEM and those from the formulas will be 15.79% in [R] and 4.54% in [L] from 1Hz to 1MHz. By decreasing the conductor radii of SC cables, the distances between the cables and the pipe will be increased, which is likely to decrease the errors. As a test case, if rB^^mm, r i =20mm, rBj=23mm, and r 4 =25mmforthe PT cable in Fig. 6.4(b) with J 2 J 3 the other parameters remaining the same, the maximum differences between the results from FEM and those from the formulas will be 13.46% in [R] and 24.18% in [L] from 1Hz to 1MHz. A similar PT cable was studied in reference[42]. The distances between the SC cables and the pipe in that case are larger than those in Fig. 6.4, but smaller than in the above test case. Accordingly, the differences between the results from FEM and those from the formulas lie also in the middle. From the above discussions, it is concluded that the formulas for the impedance calculation of PT cables give accurate results in the low frequency range. They yield reasonably accurate results in the high frequency range if the cable dimensions are small compared with the dimension of the pipe and if the cables and the pipe are not too close to each other. For the cradle arrangement, the impedances are listed in Tab.6.7 and Tab.6.8 with [^AA] [^CCl = a n ( ^ [^AB^f^BcJ* ^ke maximum differences between the results from FEM and those from the formulas are 81% in [R] and 72% in [L] from 1Hz to 1MHz. 6.3.2 The [C] calculation of P T cables [C] matrices for the two arrangements are given in Tab.6.9 and Tab.6.10, respectively. They are calculated with FEM and with the analytical formulas [23]. Both the energy method and the surface charge method discussed in Section 5.3 are used to calculate the capacitances. e,=l is used for all the regions. The meshes are similar to the ones in Fig. 6.5 except that the conductor regions are removed. Chapter 6. Case Studies in [Z] and [C] Calculations Table 6.7: [^AAK^CCD [^BBI 1*AA1 60 I^BBI PAAI 600 [^BB! I*AA1 6000 (^BB) [*AA1 60000 I^BB) t^AAl 600000 Table 6.8: #12 #22 Xll X12 L22 [%B1 FEM ana. 0.285638 0.279190 0.287547 0.279190 1.01671 0.977012 1.04002 0.977012 3.65653 3.89726 3.56767 3.89726 16.2286 16.3205 14.9878 16.3205 0.264232 0.257818 0.266144 0.257818 0.957819 0.917958 0.981136 0.917958 3.46628 3.70623 3.37744 3.70623 14.9059 15.0007 13.6651 15.0007 0.666324 0.659823 0.668125 0.659823 1.35986 1.31981 1.38306 1.31981 3.86192 4.10172 3.77298 4.10172 14.9714 15.0675 13.7305 15.0675 1.04983 1.07678 1.04277 1.07678 0.90185 0.92884 0.80478 0.92884 0.90024 0.92721 0.80317 0.92721 PEM ana. FEM ana. 0.50716 0.55370 0.40151 0.55370 0.29510 0.34546 0.27689 0.34546 0.38652 0.43298 0.37087 0.43298 0.18432 0.23464 0.16611 0.23464 0.38491 0.43137 0.36026 0.43137 0.18272 0.23304 0.16451 0.23304 PEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. 73.0357 56.6433 64.4238 56.6433 68.5463 52.1496 50.9340 52.1496 68.5544 52.1582 50.0430 52.1582 0.21712 0.25717 0.20428 0.25717 0.17656 0.22257 0.16810 0.22257 0.11126 0.15128 0.09843 0.15128 0.07407 0.12004 0.06570 0.12004 0.11053 0.15055 0.09769 0.15055 0.07407 0.12004 0.06570 0.12004 [^ABK^BC]) anc (Hz) PABI l*ACl 1*AB1 [*ACl [*AB! 6000 [*ACl [*AB! 60000 [*AC] 1*AB1 600000 ^ l^ACl °^ PT cable with a cradle arrangement L (mH/km) # (ft/km) / 600 L (mH/km) #11 (Hz) 60 °f the PT cable with a cradle arrangement # (fi/km) / t*ACl FEM ana. PEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. FEM ana. 126 #11 #12 #21 #22 £n 0.2513 0.2508 0.2446 0.2461 0.7972 0.7855 0.7149 0.7071 2.663 2.351 2.218 1.759 9.605 6.832 7.024 4.502 35.39 20.51 21.59 12.93 0.2513 0.2508 0.2446 0.2461 0.7972 0.7855 0.7149 0.7071 2.663 2.351 2.218 1.759 9.605 6.832 7.024 4.502 35.39 20.51 21.50 12.93 0.2513 0.2508 0.2446 0.2461 0.7972 0.7855 0.7149 0.7071 2.663 2.351 2.218 1.759 9.605 6.832 7.024 4.502 35.39 20.51 21.50 12.93 0.2513 0.2508 0.2446 0.2461 0.7972 0.7855 0.7149 0.7071 2.663 2.351 2.218 1.759 9.605 6.832 7.024 4.502 35.39 20.51 21.59 12.93 0.7003 0.6981 0.6115 0.5975 L\ L21 L22 0.7003 0.6981 0.6115 0.5975 0.7003 0.6981 0.6115 0.5975 0.7003 0.6981 0.6115 0.5975 0.2426 0.2429 0.1845 0.1689 0.0940 0.1046 0.0635 0.0622 0.2426 0.2429 0.1845 0.1689 0.0940 0.1046 0.0635 0.0622 0.2426 0.2429 0.1845 0.1689 0.0940 0.1046 0.0635 0.0622 0.2426 0.2429 0.1845 0.1689 0.0940 0.1046 0.0635 0.0622 0.0407 0.0651 0.0218 0.0374 0.0407 0.0651 0.0218 0.0374 0.0407 0.0651 0.0218 0.0374 0.0407 0.0651 0.0218 0.0374 0.0208 0.0534 0.0094 0.0302 0.0208 0.0534 0.0094 0.0302 0.0208 0.0534 0.0094 0.0302 0.0208 0.0534 0.0094 0.0302 2 Chapter 6. Case Studies in [Z] and [C] Calculations 127 Table 6.9: [C] of the PT cable with the triangle arrangement (/zF/km) [C] found with FEM using energy method 0.109652 -0.109652 0 0 0 0 -0.109652 0 0.243091 0 0 0.109652 -0.040218 -0.109652 0 0 -0.040218 0 0 -0.040218 -0.109652 0.342047 0 -0.014692 0 0 0 0 0.109652 -0.109652 0 -0.040218 0 -0.014692 -0.109652 0.342052 [C] found with FEM using surface charge method 0.109442 -0.109592 0 0 0 0 -0.109442 0 0.243030 0 0 0.109442 -0.040217 -0.109593 0 0 -0.040217 0 0 -0.040217 -0.109442 0.341983 0 -0.014692 0 0 0 0 0.109442 -0.109593 0 -0.040217 0 -0.014692 -0.109442 0.341988 [C] found with analytical formulas 0.109795 -0.109795 0 0 0 0 -0.109795 0 0.209431 0 0 0.109795 -0.039429 -0.109795 0 0 -0.039429 0 0 -0.039429 -0.109795 0.243418 0 -0.022992 0 0 0 0 0.109795 -0.109795 0 -0.039429 0 -0.022992 -0.109795 0.243418 Table 6.10: [C] of the PT cable with the cradle arrangement (/iF/km) [C] found with FEM using energy method 0.109653 -0.109653 0 0 0 0 -0.109653 0 0.348885 0 0 0.109651 -0.053153 -0.109651 0 0 -0.007603 0 0 -0.053153 -0.109651 0.386509 0 -0.053166 0 0 0 0 0.109653 -0.109653 0 -0.007603 0 -0.053166 -0.109653 0.348895 [C] found with FEM using surface charge method 0.109442 -0.109593 0 0 0 0 -0.109442 0 0.348820 0 0 0.109440 -0.053152 -0.109590 0 0 -0.007603 0 0 -0.053152 -0.109440 0.386443 0 -0.053165 0 0 0 0 0.109442 -0.109593 0 -0.007603 0 -0.053165 -0.109442 0.348831 [C] found with analytical formulas 0.109795 -0.109795 0 0 0 0 -0.109795 0 0.243352 0 0 0.109795 -0.057283 -0.109795 0 0 -0.007962 0 0 -0.057283 -0.109795 0.267446 0 -0.057283 0 0 0 0 0.109795 -0.109795 0 -0.007962 0 -0.057283 -0.109795 0.243352 Chapter 6. Case Studies in [Z] and [C] Calculations 6.4 [Z] and [C] Calculations of Sector-Shaped Cables 6.4.1 The [Z] calculation of sector-shaped cables 128 There is no analytical formula for the impedance calculation of sector-shaped cables of the type on the right in Fig. 1.1(a). The sheath indicated in that figure may or may not exist. For a sector-shaped cable with a sheath, Ametani suggested that its impedances can be calculated approximately with the formulas for PT cables[40]. If Lc and Sc represent the contour length and cross-section area of a sector-shaped core conductor in Fig. 1.1(a), respectively, the conductor can be transformed into an equivalent circular conductor with its outer and inner radii as [40] outer radius inner radius r = TJ— 27T B = )jr — B (6-2) (6-3) Reference [40] did not mention, however, how to determine the position of the equivalent conductor inside the sheath. The impedances of the sector-shaped cable in Fig. 6.8(a) are calculated with FEM and with the analytical formulas for PT cables using the equivalent radii. For one sector-shaped conductor, 5'c=300mm and Lc=70.834mm, with the equivalent radii be2 ing r4=5.622mm and rB=11.274mm. It is assumed that the three equivalent circular j conductors are 13.1mm away from the center of the cable. The sheath is used as the return path and the earth is not considered. A=Q is assumed on the boundary located at r in Fig. 6.8(a). The FEM mesh around the cable at 60kHz is plotted in Fig. 6.8(b). 4 The impedances of the cable are listed in Tab.6.11. Due to symmetry, Zi 1 = ^ 2 2 = ^ 3 3 and ^ 1 2 = ^ 1 3 = ^ 2 3 . "app" in the table represents the results found with the approximate formulas for PT cables using equivalent radii. For a loop current 1+jO A between the upper conductor and the sheath, with the upper conductor current being 1 A, the |A| Chapter 6. Case Studies in [Z] and [C] Calculations 129 t =4.255 mm (a) cable data (b) FEM mesh around cable at 60kHz Figure 6.8: A sector-shaped cable distributions at 60Hz and 60kHz are given in Fig. 6.9(a) and (b), respectively. The plotting area is 100mm by 100mm. Table 6.11: [Z] of the sector-shaped cable R (Jtykm) / (Hz) 6 60 600 6000 60000 600000 FEM app FEM app FEM app FEM app FEM app FEM app Rn R12 2.84006 2.83996 2.84987 2.84334 2.96756 2.93723 3.52505 3.52328 5.15051 5.89100 15.9166 18.8742 2.78246 2.78239 2.78317 2.78105 2.78308 2.75504 2.69462 2.58640 2.88824 2.22645 8.82314 6.42831 L (/iH/km) 232.020 196.081 220.673 191.259 156.797 164.736 120.400 125.467 105.090 108.540 99.0031 100.191 40.4454 25.0563 40.1914 27.3958 35.7225 36.4771 38.3245 47.2753 40.6921 52.1332 38.0609 51.1033 The results show that in the low frequency range the differences in [R] between FEM and the formulas are very small because the equivalent circular conductor has the same cross-section area as a sector-shaped conductor. The differences in [£], however, are very large at low frequencies because the magneticfieldis not confined within the cable. This is clearly indicated by the \A\ distribution at 60Hz in Fig. 6.9(a). In the high frequency range, large differences in [R] and in Lu can be observed, while the differences in L n Chapter 6. Case Studies in [Z] and [C] Calculations 130 |4| = 1.8777 + 14.298- (i-l) \At\ = 0.04552 + 5.7462- (i-l) (a) \A\ at 60Hz (/«V/km) (b) \A\ at 60kHz (/isV/km) Figure 6.9: \A\ distribution in the sector-shaped cable caused by the loop current become very small. The magnetic field becomes confined within the cable as shown by the \A\ distribution at 60kHz in Fig. 6.9(b). The maximum differences are 41.25% in [R] and 62.11% in [L] from 1Hz to 1MHz. If the sheath is made of steel with high permeability the magnetic field will be confined within the cable even at low frequencies, as seen in Fig. 6.6(a) in the preceding section. The formulas may give reasonably accurate results for [L] in the low frequency range as well. If /z,=500 for the sheath, the difference in [L] between FEM and the formulas is less than 5% if the frequency is less than 2kHz, and the maximum differences are 23.86% in [R] and 17.57% in [L] from 1Hz to 1MHz. 6.4.2 The [C] calculation of sector-shaped cables Capacitances of the cable are calculated with FEM and with the formula for PT cables. The results are listed in Tab.6.12 for e =l. If the potential of the upper conductor is IV r and the other conductors and the sheath have zero potentials, the potential distribution Chapter 6. Case Studies in [Z] and [C] Calculations 131 Table 6.12: [C] of the sector-shaped cable (pF/km) method C l 1(^22 ^ 3 3 ) Cl2(Cl3>(723) FEM (energy method) FEM (surface charge method) analytical formula 0.147049 0.147323 0.181528 -0.0401938 -0.0402814 -0.0620687 found with FEM is plotted in Fig. 6.10. Fringe effects around the corners of the upper conductor can be observed. Because of numerical round-off errors, 4>% given in Fig. 6.10 starts from -0.000003 instead of zero. <t>i = -0.000003 + 0.05- (i-l) V Figure 6.10: Potential distribution in the sector-shaped cable 6.5 The Calculation of Internal Resistance of Stranded Conductors Stranded conductors are used for overhead transmission lines. In such lines, the external inductance is much larger than the internal inductance of the conductors. The internal inductance is therefore ignored, or calculated approximately. However, the internal resistance of the conductors is of great importance. By ignoring the spiralling effect, the Chapter 6. Case Studies in [Z] and [C] Calculations 132 strands of a stranded conductor can be assumed as parallel circular subconductors. With such an assumption, analytical formulas can be derived for the calculation of the internal resistance of stranded conductors. Galloway, Shorrocks, and Wedepohl derived the following formula in 1964 for the internal resistance in the high frequency range[4] where n is the number of the outer strands, r is the radius of the outer strands, 6 is the real penetration depth defined in (3.45), a is the conductivity of the stranded conductor, and Kj=2.25 is a factor found from measurements in an electrolytic tank. The formula shall be called the GSW formula for simplicity. Borges da Silva suggested in 1979 that a stranded conductor could be replaced with a circular conductor having an equivalent radius[21]. The internal resistance can then be calculated with the formula for the impedance calculation of circular conductors, using the equivalent radius given as r ea e * = r.(0.92122679 1 V 4.3856939n - 2.3071869n -1.2479854^ e 2 (6.5) K } where r is the outer radius of the stranded conductor, n is the same as in (6.4). This e method shall be called Borges da Silva's formula. Arizon et al applied the conductor subdivision method to calculate the internal impedances of stranded conductors in 1987[38]. In that paper differences were reported between the results from the subdivision method and those from the GSW formula, which become small if lfy=1.6 is used in the formula instead of ify=2.25. With the assumption that the strands are parallel, FEM can be applied to the calculation of internal resistance of stranded conductors. A two-layer stranded conductor is given in Fig. 6.11(a). It is assumed that the conductor is surrounded by the air. A circular boundary away from the conductor is set up, with A=0 on the boundary. The Chapter 6. 133 Case Studies in [Z] and [C] Calculations real part of the impedance for such a single conductor system would be the internal resistance of the conductor. Due to symmetry, only one twelfth of the conductor is used in the FEM solution, as shown in Fig. 6.11(b). The FEM mesh around the conductor at 60kHz is plotted in Fig. 6.11(c). (a) data (b) solution region (c) FEM mesh at 60kHz (partial) Figure 6.11: A two-layer stranded conductor The internal resistance of the conductor in Fig. 6.11 is given in Tab.6.13. The maximum differences between the results from FEM and those from the other methods are Table 6.13: The internal resistance of the two-layer stranded conductor / (kHz) FEM 43 1.0446 80 1.4160 100 1.5905 1.7951 130 max diff with FEM Rc (A/km) GSW formula subdivision #,=2.25 #,=1.6 method 1.4524 1.0328 1.3479 1.9811 1.4088 1.6305 2.2149 1.5751 1.8187 2.5254 1.7959 2.1126 29.04% 40.68% 1.13% Borges da Silva's formula 1.0215 1.3864 1.5479 1.7622 2.68% given at the bottom of Tab.6.13. With r =13mm and n=12, r ,=11.7171mm is found e e from (6.5). It can be seen that differences between FEM and Borges da Silva's formula are small. There are large differences between FEM and the GSW formula with Kf=2.25. If Kf=1.6 is used in the GSW formula, the differences between FEM and the formula become insignificant, as shown in Tab.6.13. Chapter 6. Case Studies in [Z] and [C] Calculations 134 If the current in the stranded conductor is 1 A, the \J\ distribution at 60Hz and 60kHz are plotted in Fig. 6.12. At 60Hz the currents are almost evenly distributed over \Ji\ = 2455.8 + 8.21- (i-l) \J \ = 3520.7- (i-l) (a) |J| at 60Hz (A/m ) (b) \J\ at 60kHz (A/m ) t 2 2 Figure 6.12: \J\ distribution in the stranded conductor the conductor, while at 60kHz a strong skin effect can be observed. If the stranded conductor in Fig. 6.11(a) has only one layer, the corresponding internal resistance is given in Tab.6.14, for n=6, r =7.8mm, and r ,=6.8578mm. The differences e e Table 6.14: The internal resistance of the one-layer stranded conductor / (kHz) FEM 20 1.2436 43 1.7882 80 2.4094 100 2.7013 130 3.0466 max diff with FEM Rc (n/km) GSW formula subdivision K =2.25 #,=1.6 method 1.7335 1.2327 1.7122 2.5418 1.8075 1.8466 3.4669 2.4654 2.4815 3.8761 2.7564 2.8160 4.4195 3.1427 3.2965 37.68% 45.06% 3.15% f Borges da Silva's formula 1.2246 1.7688 2.3920 2.6677 3.0338 1.53% are very small among the results of FEM, Borges da Silva's formula, and the GSW formula with Kf=1.6. Good agreement is observed between FEM and the conductor subdivision method for the listed frequencies, except at 20kHz. For three-layer and four-layer stranded conductors with the same D as in Fig. 6.11(a), n will be 18 and 24, respectively, and r will be 18.2mm and 23.4mm, respectively, with r = 16.5286mm e eq Chapter 6. 135 Case Studies in [Z] and [C] Calculations for the three-layer stranded conductor and r = 21.3293mm for the four-layer stranded eq conductor. The internal resistance is listed in Tab.6.15 and Tab.6.16, respectively. Table 6.15: The internal resistance of the three-layer stranded conductor / (kHz) FEM 20 0.50658 60 0.87271 1.1174 100 140 1.3169 180 1.4950 220 1.6610 max diff with FEM Rc (O/km) GSW formula #,=2.25 0.69338 1.2010 1.5505 1.8345 2.0802 2.2997 39.30% #,=1.6 0.49307 0.85403 1.1025 1.3045 1.4792 1.6353 2.67% Borges da Silva's formula 0.49425 0.84906 1.0934 1.2920 1.4637 1.6172 2.71% Table 6.16: The internal resistance of the four-layer stranded conductor / (kHz) FEM 20 0.39099 60 0.67642 100 0.86673 140 1.0229 180 1.1627 220 1.2934 max diff with FEM Rc (n/km) GSW formula tf/=2.25 If) =1.6 0.53337 0.37929 0.92383 0.65694 1.1927 0.84811 1.4112 1.0035 1.1379 1.6001 1.7690 1.2580 2.99% 37.95% Borges da Silva's formula 0.38133 0.65630 0.84563 0.99951 1.1326 1.2515 3.24% It should be noted that even though the GSW formula with a factor of 2.25 does not agree with the other methods, it does seem to come closer tofieldtests[20]. The reason may be that the factor 2.25 takes the spiralling effect into account. 6.6 Summary In this chapter FEM is applied to the parameter calculation of buried or tunnel installed multiphase cables, of three-phase PT cables and sector-shaped cables, and of stranded Chapter 6. Case Studies in [Z] and [CJ Calculations 136 conductors. The eaxth region reduction technique discussed in Chapter 4 is also used in the impedance calculation of buried or tunnel installed.cables. The results show that for typical earth resistivities Pollaczek's formula gives reasonably accurate self and mutual impedances for multiphase buried cables, and for tunnel installed cables using approximate inner earth radii. Good agreement between the proposed technique and the conventional FEM are shown by the results. For PT cables, the analytical formulas may introduce errors at high frequencies because the influence of non-current carrying conductors on the fields inside the pipe and within the pipe conductor are not considered in the formulas. For sector-shaped cables comparisons show that reasonably accurate self and mutual resistances can be obtained in the low frequency range with the formulas for PT cables by using equivalent core conductor radii suggested by Ametani. Accurate self and mutual inductances cannot be obtained in the low frequency range with the formulas if the sheath is non-magnetic or if there is no sheath. In the high frequency range, large differences between the results from FEM and those from the formulas for both resistances and inductances are observed. For the internal resistance of stranded conductors, good agreement is obtained among the results from FEM, Borges da Silva's formula, and the GSW formula with Kf=1.6. There is also Good agreement between the results from FEM and those from the conductor subdivision method for a one-layer stranded conductor at />43kHz. With respect to the factor in the GSW formula, good agreement was reported between the GSW formula with ify=2.25 and field tests[20]. This may be due to the fact that the factor 2.25 takes the spiralling effect into account. Chapter 7 Conclusions and Recommendations for Future Work In this thesis, thefiniteelement method (FEM) was applied to the parameter calculation of underground power cables. The parameters most often needed in power system analysis are the series impedances and shunt capacitances. The principal equations describing the quasi-magnetic fields and static electricfieldswere solved with FEM based on the Galerkin technique. Quadratic isoparametric elements as well as high-order simplex elements were studied. A technique based on the perturbation concept was proposed to reduce the solution region in the earth. The parameters of shallowly buried or tunnel installed multiphase single core (SC) cables, pipe-type (PT) cables, sector-shaped cables, and stranded conductors were calculated with FEM. The major conclusions are: 1. The Js method and the loss-energy method derived in the thesis for the [Z] calculation from thefieldsolution gave the same results. The loss-energy method is time-consuming and requires a complete field solution. Much less computation is needed with the Js method, and a completefieldsolution is not required. 2. Quadratic isoparametric elements proved to be more efficient than high-order simplex elements, in both the impedance calculation and the capacitance calculation. For the same error tolerance, isoparametric elements can have larger span angles. The low accuracy of simplex elements for large span angles cannot be improved by increasing the order of polynomials of the elements. Fewer elements are needed in meshing a circular region with isoparametric elements. 137 Chapter 7. Conclusions and Recommendations for Future Work 138 3. Accurate impedances could be obtained for deeply buried SC cables with the conventional FEM if the field truncation boundary is at least ZS away from the cables e {n > 36e), and if the earth is divided with the pattern 10 , 10 s, 10 3, 10 n n+ N+ n+1 in each decade in the radial direction. For shallowly buried SC cables, the iterative results showed that rj > 12S is required to make the differences in [Z] from two e consecutive steps less than 0.5%. 4. Good agreement was obtained between the proposed technique and the conventional FEM for shallowly buried SC cables. As shown in Section 4.5.2, for rj,=5m, with />=100f2m and with r varying between 24mm and lm, the maximum errors with e e the proposed technique are less than 8% in [R] and less than 2% in [L]. The earth solution region is reduced significantly with the proposed technique in the low frequency range, and CPU time is saved if partial earth return currents required for the technique are calculated only once. 5. Comparisons between Pollaczek's formula and the FEM for a single-phase shallowly buried SC cable in Chapter 4 showed that accurate results were obtained with Pollaczek's formula when r /S is small. For r =24mm, the maximum differe e e ences between the results from Pollaczek's formula and from the FEM are less than 1% in [R] if r /6 <0.03 and in [L] if r /8 <0.095. For large r /«5 , the differences e e e e e e become large. With typical ranges of p and r , however, the maximum differences e e between the two approaches from 1Hz to 1MHz are reasonably small. With p e varying between lOOOfim to l£2m and with r varying between 24mm to 250mm, e the maximum differences are less than 21% in [R] and less than 9% in [L]. Results also showed that Pollaczek's formula could be applied to find the impedances of tunnel installed cables as well, if an approximate r is used. For typical earth resise tivities, Pollaczek's formula gave reasonably accurate self and mutual impedances Chapter 7. Conclusions and Recommendations for Future Work 139 for multiphase buried cables, and for tunnel installed cables with an approximate r,. 6. For a PT cable, the magneticfluxdistribution in the pipe and the current densitydistribution within the pipe conductor are significantly influenced by the presence of non-current carrying conductors at high frequencies. This influence is ignored in approximate formulas, which therefore produce noticeable errors at high frequencies. 7. For a sector-shaped cable with a magnetic sheath, the the analytical formulas suggested by Ametani produce reasonably accurate resistances and inductances in the low frequency range. If the sheath is non-magnetic or nonexistent, then the inductance has a large error. In the high frequency range, these formulas are too inaccurate. 8. For the calculation of the internal resistance of stranded conductors, the spiralling effect was ignored by assuming that all strands are in parallel. Close agreement was obtained among the results from FEM, Borges da Silva's formula, and the GSW formula with ify=1.6, for stranded conductors with one to four layers. It has been reported, however, that the GSW formula with Kf=2.25 comes closer to field tests. This may be due to the fact that the factor 2.25 takes the spiralling effect into account. There was also good agreement between the results from FEM and those from the conductor subdivision method for a one-layer stranded conductor at />43kHz. 9. For the capacitance calculation from the field solution, the energy method is simple and easy to implement. It has less stringent requirements on the mesh and has higher accuracy than the surface charge method. The surface charge method, Chapter 7. Conclusions and Recommendations for Future Work 140 however, is faster, and its accuracy can be improved with a finer mesh near the conductor surfaces. Future research could be conducted in the following areas: 1. The mesh should be automatically generated for cable geometries, for user-friendly interfaces with the FEM. The auto-mesh program developed for this project is not flexible enough, and is incomplete. 2. More field tests are needed to compare the calculated impedances against measured impedances. 3. If the conductance of insulating materials is available from tests and is to be considered as well, then the FEM should be modified to include it. 4. the FEM program could be used as a verification tool for developing simpler approximate formulas for PT cables, sector-shaped cables, and other types of cables. References [1] J. R. Carson, "Wave Propagation in Overhead Wires, with Ground Return," Bell Syst. Tech. Jour., vol.5, 1926, pp.539-554. [2] F. Pollaczek, "Uber das Feld einer unendlich langen wechselstromdurchfiossenen Einfachleitung," E.N.T., 1926, Band 3 (Heft 9), pp.339-360. [3] S.A. Schelkunoff,"The Electromagnetic Theory of Coaxial Transmission Line and Cylindrical Shells," Bell System Tech. Jour., vol.13, 1934, pp.522-579. [4] R. H. Galloway, W. B. Shorrocks, and L. M. Wedepohl, "Calculation of Electrical Parameters for Short and Long Polyphase Transmission Lines," Proc. IEE, vol.111, Dec. 1964, pp.2051-2059. [5] A. H. Stroud and Don Secrest, Gaussian Quadrature Formulas. Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1966. [6] L. M. Wedepohl and S. E. T. Mohamed, "Multiconductor Transmission Lines — Theory of Natural Modes and Fourier Integral Applied to Transient Analysis," Proc. IEE, vol.116, Sept. 1969, pp.1553-1563. [7] P. Silvester, "High-Order Polynomial Triangular Finite Elements for Potential Problems," Int. J. Engng Sci., vol.7, 1969, pp.849-861. [8] J. A. Tegopoulos and E. E. Kriezis, "Eddy Current Distribution in Cylindrical Shells of Infinite Length due to Axial Currents - Part II: Shells of Finite Thickness," IEEE Trans, on PAS, vol.PAS-90, May 1971, pp.1287-1294. [9] E. Comellini, A. Invernizzi, and G. Manzoni, "A Computer Program for Determining Electrical Resistance and Reactance of any Transmission Line," IEEE Trans. PAS, vol.PAS-92, Jan./Feb. 1973, pp.308-314. 141 References 142 [10] L. M. Wedepohl and D. J. Wilcox, "Transient Analysis of Underground PowerTransmission Systems — System Model and Wave-Propagation Characteristics," Proc. IEE, vol.120, Feb. 1973, pp.253-260. [11] H. C. Martin and G. F. Carey, Introduction to Finite Element Analysis. McGraw- Hill Book Company, 1973. [12] M. V. K. Chari, "Finite Element Solution of the Eddy-Current Problem in Magnetic Structures," IEEE Trans. PAS, vol.PAS-93, Jan. 1974, pp.62-72. [13] G. W. Brown and R. G. Rocamora, "Surge Propagation in Three-Phase Pipe-Type Cables. Part I - Unsaturated Pipe," IEEE Trans. PAS, vol. PAS-95, Jan./Feb. 1976, pp.89-95. [14] M. V. K. Chari and Z. J. Csendes, "Finite Element Analysis of the Skin Effect in Current Carrying Conductors," IEEE Trans, on Magnetics, vol.MAG-13, Sept. 1977, p.1125-1127. [15] A. Jennings, Matrix Computation for Engineers and Scientists. John Wiley & Sons Ltd., 1977. [16] E. M. Deeley and E. E. Okon, " An Integral Method for Computing the Inductance and the A.C. Resistance of Parallel Conductors," Intl. J. Numerical Method Eng., vol.12, 1978, pp.625-634. [17] R. Lucas and S. Talukdar, "Advances in Finite Element Techniques for Calculating Cable Resistances and Inductances," IEEE Trans. PAS, vol.PAS-97, May/June 1978, pp.875-883. [18] L. M. Wedepohl and A. E. Efthymiadis, "Wave Propagation in Transmission Lines over Lossy Ground: a New, Complete Field Solution," Proc. IEE, Vol.125, June 1978, pp.505-510. [19] A. E. Efthymiadis and L. M. Wedepohl, "Propagation Characteristics of InfinitelyLong Single-Conductor Lines by the Complete Field Solution Method," Proc. IEE, 143 References Vol.125, June 1978, pp.511-517. [20] S. Cristina and M. D'Amore, "Propagation on Polyphase Lossy Lines: A New Parameter Sensitivity Model," IEEE Trans. PAS, vol.PAS-98, Jan./Feb. 1979, pp.3544. [21] J. F. Borges da Silva, "The Electrostatic Field Problem of Stranded and Bundle Conductors Solved by the Multipole Method," Electricidade, No.142, Mar./Apr. 1979, ppl-11. [22] W. T. Weeks, L. L. Wu, M. F. McAllister, and A. Singh, "Resistive and Inductive Skin Effect in Rectangular Conductors," IBM J. RES. DEVELOP, vol.23, Nov. 1979, p 652-660. P [23] A. Ametani, "A General Formulation of Impedance of Cables," IEEE Trans. PAS, vol.PAS-99, May/June 1980, pp.902-909. [24] J. M. Schneider and S. J. Salon, "A Boundary Integral Formulation of the Eddy Current Problem," IEEE Trans. MAG, vol.MAG-16, Sept. 1980, pp.1086-1088. [25] W. D. Humpage, K. P. Wong, T. T. Nguyen, and D. Sutanto, "z-Transform Electromagnetic Transient Analysis in Power Systems," Proc. IEE-C, Vol.127, Nov. 1980, pp.370-378. [26] Al Konrad, "The Numerical Solution of Steady-State Skin Effect Problems — An Integrodifferential Approach," IEEE Trans. MAG, vol.MAG-17, Jan. 1981, pp.1148-1152. [27] E. L. Wachspress, "High-Order Curved Finite Elements," Intl. J. Numerical Method Eng., vol.17, 1981, pp.735-745. [28] A. Konrad, "Integrodifferential Finite Element Formulation of Two Dimensional Steady-State Skin Effect Problems," IEEE Trans. MAG, vol.MAG-18, Jan. 1982, pp.284-292. References 144 [29] J. R. Marti, "Accurate Modelling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulations," IEEE Trans. PAS, vol.PAS-101, Jan. 1982, pp.147-157. [30] S. J. Salon and J. M. Schneider, "A Hybrid Finite Element-Boundary Integral Formulation of the Eddy Current Problem," IEEE Trans. MAG, vol.MAG-18, Mar. 1982, pp.461-466. [31] J. Weiss and Z. J. Csendes, "A One-Step Finite Element Method for Multiconductor Skin Effect Problems," IEEE Trans, on PAS, vol.PAS-101, Oct. 1982, pp.3796-3803. [32] J. Weiss, V. K. Garg, and E. Sternheim, "Eddy Current Loss Calculation in Multiconductor Systems," IEEE Trans. MAG, vol.MAG-19, Sept. 1983, pp.2207-2209. [33] W. M. Rucker and K. R. Richter, "Calculation of Two-Dimensional Eddy Current Problems with the Boundary Element Method," IEEE Trans. MAG, vol.MAG-19, Nov. 1983, pp.2429-2432. [34] P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers. Cam- bridge University Press, 1983. [35] S. J. Salon, "The Hybrid Finite Element-Boundary Element Method in Electromagnetics," IEEE Trans. MAG, vol.MAG-21, Sep. 1985, pp.1829-1834. [36] J. Poltz and E. Kuffel, "Application of Boundary Element Techniques for 2D EddyCurrent Problems," IEEE Trans. MAG, vol.MAG-21, Nov. 1985, pp.2254-2256. [37] P. E. Allaire, Basics of the Finite Element Method — Solid Mechanics, Heat Trans- fer, and Fluid Mechanics. Wm. C. Brown Publishers, Dubuque, Iowa, 1985. [38] P. de Arizon and H. W. Dommel, "Computation of Cable Impedances Based on Subdivision of Conductors," IEEE Trans, on Power Delivery, vol.PWRD-2, Jan. 1987, pp.21-27. References 145 [39] L. Marti, "Simulation of Transients in Underground Cables with FrequencyDependent Modal Transformation Matrices," IEEE Trans, on Power Delivery, vol.PWRD-3, July 1988, pp.1099-1110. [40] A. Ametani, "Further Improvements of CABLE CONSTANTS and An Investigation of Cable Problems," EMTP NEWS, vol.1, No.4, Dec. 1988, pp.4-14. [41] S. Cristina and M. Feliziani, "A Finite Element Technique for Multiconductor Cable Parameters Calculation," IEEE Trans, on Magnetics, vol.MAG-25, July 1989, pp.2986-2988. [42] Y. Yin and H. W. Dommel, "Calculation of Frequency-Dependent Impedances of Underground Power Cables with Finite Element Method," IEEE Trans, on Magnetics, vol.MAG-25, July 1989, pp.3025-3027. [43] M. Rioual, "Measurements and Computer Simulation of Fast Transients through Indoor and Outdoor Substations," IEEE Trans, on Power Delivery, vol.PWRD-5, Jan. 1990, pp.117-123. [44] A. Konrad and P. Silvester, "Triangular Finite Elements for the Generalized Bessel Equation of Order m," Int. J. Numer. Meth. Eng., vol.7, 1973, pp.43-55. [45] B. Engquist, A. Greenbaum, and W. D. Murphy, "Global Boundary Conditions and Fast Helmholtz Solvers," IEEE Trans, on Magnetics, vol.MAG-25, July 1989, pp.2804-2806. Appendix A Integral matrices [Q^] and [Ts] of simplex elements Table A.l: [Q^] and [Ts] of the 1st order simplex element (a) (b) [Ts] common denominator = 2 common denominator = 12 0 symm. 0 1 0-11 2 symm. 12 1 1 2 Table A.2: [QW] and [T ] of the 2nd order simplex element s (a) (b) [T ] S common denominator = 6 0 0 8 0 -8 0 0 0 0 0 0 common denominator = 180 6 0 32 symm. 0 16 32 -1 0 -4 6 -4 16 16 0 32 -1 -4 0 -1 0 6 symm. 8 0 3 0 -4 8 0 1 -4 3 146 Appendix A. Integral matrices [QW] and [Ts] of simplex elements Table A.3: [QM] and [T ] of the 3rd order simplex element s (a) (b) [Ts] common denominator = 6720 [Q ] (1) common denominator = 80 0 0 135 symm. 0 -135 135 0 -27 27 135 0 0 0 -162 324 0 27 -27 27 -162 135 0 3 -3 3 0 -3 34 0 0 0 0 0 0 -54 135 0 0 0 0 0 0 27 -108 135 0 -3 3 -3 0 3 -7 27 -54 76 18 540 18 270 0 -189 36 162 0 -135 11 0 27 -135 27 -54 11 27 34 symm. 540 -135 162 -189 27 -54 -135 0 540 162 1944 -54 162 540 18 36 27 270 162 -135 -135 162 270 27 36 18 76 18 540 0 -189 11 0 Table A.4: [QW] and [T ] of the 4th order simplex element s (a) [<2 ] (common denominator = 1890) (1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3968 -3968 -1440 0 1440 640 0 0 -640 -80 0 0 0 80 3968 1440 0 -1440 -640 0 0 640 80 0 0 0 -80 4632 -5376 10752 744 -5376 4632 -1248 1536 -288 768 -1536 768 768 -1536 768 -288 1536 -1248 80 -160 80 -128 256 -128 96 -192 96 -128 256 -128 80 -160 80 symm, 3456 -4608 10752 1536 -7680 10752 -384 1536 -4608 240 -160 -160 -256 256 256 192 -192 -192 -256 256 256 80 -160 -160 3456 80 -256 192 -256 240 705 -1232 884 -464 107 3456 -3680 5592 1920 -3680 -464 884 3456 -1232 705 2560 160 290 (b) [T ] (common denominator = 56700) s 290 160 2560 160 1280 -80 -1280 160 1280 -80 -960 0 768 -160 256 -160 -256 0 512 -27 0 -112 512 -12 64 -112 256 -27 -112 2560 -960 1280 -1280 512 -256 256 768 -112 256 64 512 0 3168 384 10752 48 384 -1280 256 384 -1536 -768 -1536 64 256 -80 -160 -960 -256 48 -768 64 -256 -12 -160 symm. 3168 64 -768 384 -1280 -12 64 48 -960 -80 2560 1280 10752 -256 -1536 10752 256 -256 1280 160 160 -160 1280 1280 256 -960 384 384 512 256 1280 -112 -160 160 2560 -112 512 -960 1280 160 290 160 -80 0 -27 2560 -1280 768 0 3168 -1280 -80 (a) (common denominator = 290304) 0 743750 0 -743750 743750 0 -405000 405000 1072500 0 0-1150000 2300000 0 0 405000 -405000 77500-1150000 1072500 0 260000 -260000 -523750 587500 -63750 892500 0 0 0 293750 -587500 293750-1106250 2512500 0 0 0 293750 -587500 293750 300000-1706250 2512500 0 -260000 260000 -63750 587500 -523750 -86250 300000-1106250 892500 0 -98750 98750 237500 -287500 50000 -263125 346875 -121875 38125 0 0 0 -150000 300000 -150000 212500 -525000 412500 -100000 0 0 0 12500 -25000 12500 112500 -112500 -112500 112500 0 0 -150000 300000 -150000 -100000 412500 -525000 212500 0 0 98750 -98750 50000 -287500 237500 38125 -121875 346875 -263125 0 11850 -11850 -20225 31250 -11025 -4600 -15625 31250 -11025 0 0 0 20625 •41250 20625 -14375 61250 -79375 32500 0 0 -5000 10000 -5000 30000 -92500 95000 -32500 0 0 0 0 -5000 10000 -5000 -325O0 95000 -92500 30000 0 0 20625 -41250 20625 32500 -79375 61250 -14375 0 0 -11850 11850 -11025 31250 -20225 -11025 31250 -15625 -4600 577500 -818750 2050000 393750-1837500 2887500 -193750 800000-1837500 2050000 41250 -193750 393750 -818750 577500 58725 -62500 -15625 31250 -11850 -83125 146250 14375 -135000 57500 83125 -177500 1250 197500 -104375 -104375 197500 1250 -177500 83125 575O0 -135000 14375 146250 -83125 -11850 31250 -15625 -62500 58725 (b) [771 (commoi 19160064) 53244 45270 563500 45270 281750 563500 -36720 -367250 -281750 34200 392500 392500 -36720 -281750 -367250 17730 295250 220750 -37350 3750 -146250 -37350 -146250 3750 17730 220750 295250 -880 -152125 -95125 33700 -107500 30000 24900 -52500 -52500 33700 30000 -107500 -880 -95125 -152125 4747 -880 13935 13935 -95125 -15750 7940 -1750 -55500 7940 -55500 -1750 13935 -15750 -95125 4747 13935 -880 846000 -55000 2700000 113000 -55000 846000 -516500 195000 -52000 202500 -600000 -210000 -210000 -600000 202500 -52000 195000 -516500 295250 -107500 -1750 195000 450000 120000 -22500 150000 -22500 120000 450000 195000 -1750 -107500 295250 17730 33700 7940 220750 30000 -55500 -52000 120000 108000 108000 120000 -52000 -55500 30000 220750 7940 33700 17730 0 — In* symm. > Cn 55 99402 -188125 172500 -130625 58750 -11902 577500 -745625 1282500 570000-1148750 1282500 -272500 570000 -745625 577500 58750 -130625 172500 -188125 99402 fcr* a Cn XT O H 846000 202500 2925000 -22500 -450000 108000 -22500 -367250 3750 -55000 -600000 •210000 -450000 120000 150000 -55500 -52500 -36720 -37350 -281750 -146250 113000 -210000 -52000 -22500 -1750 -52500 7940 24900 &* 2925000 202500 846000 -52500 -55500 150000 120000 -450000 -210000 -600000 -55000 3750 •367250 24900 7940 -52500 -1750 -22500 -52000 -210000 113000 -146250 -281750 -37350 -36720 n B *—• symm. 392500 -146250 30000 -15750 45270 281750 -281750 220750 -95125 13935 2700000 -600000 2925000 450000 -600000 2700000 30000 -146250 392500 563500 34200 -37350 33700 13935 392500 3750 -107500 -95125 -55000 202500 195000 220750 195000 202500 -55000 -281750 -107500 3750 392500 281750 33700 -37350 34200 45270 n X £L n B 53244 45270 -36720 17730 -880 4747 n 563500 -367250 846000 295250 -516500 846000 •152125 295250 •367250 563500 -880 17730 -36720 45270 P 53244 Appendix A. Integral matrices [Q^] and jTs] of simplex elements 149 Table A.6: [QW] and [T ] of the 6th order simplex element s S§ s III ^9 I 888 ill S Si ill? ll » ?n z s lis M ill & E li IH «o»« saa :2SG ' *9 T till lis S J 3 H S SS - 8 : go r> i 5 9 JS8 ON i I-i O «-» c3 C 3 J S S I sa 1 Of W N n o o evoo rt oe s?s « o3oo8r« so e as S8 X 8 s o c <o S S 88 S 5! ' R S a§ 2 S 3 3 S cn* *f 3 o o o 8 8 3 r S3 5 8 sa*s I8 2 3 S P o3 SS 8 i lap i §s list 8 8 S 0 O O a R oooooooooooo ©OO' T S 3 T i i £§ i i 11 J I I c = S S 11 11 I ' M V R a a' - ^ T s a p «i 1 p i P C S * : 5 - §i m m | i s 11 § 8 IS8 §8 •* £ J£ o o o o si a OOO OOOOI O i1 o o oo 11 s i | 8 11 s »* 11$ ? s"? aS 5 8 a s aP S ! a a a" l i a § sI§ 111 8 5 1 S "7 S3 1 I li 15 n ?i § S 80 X lis <»"i ft yj I «o § s i • ° ° a R s§s a ill 88 S1S IIS s a' 5 "? S s m o « «o Ro o o e £5 0 o 2 SSSR- R ss s°°°a 0 a T T sspss^sa^aac? » * 8 s te ^ a • 2 a Tm tn oo & §3 a s s 223 Ills 1 s s sa s 8 S§1 33 a s; a s a §s : S saaj si?** 151 S s s |a i I i 1111 g s 8 I I g 9 H I 8 ag« : s a 8 B 5 SI S :B P il = § § 1 5tt« 5 8 a g1 1 1 1 ! s s* s 5g T g 8 S 8 f s a s sp a a r« <n a «• s s : !s is g I 8 | x = e a a 2 a: § s| ISO * l a s o so r 5 3 8 X saR '•as aa2 s p a a ss x a S P £ !8 I a S R S S S ! 8 tsSI R 5 3 $ a s a s s a" 2 S S S S 2 S « a p; lis 3 a o c •3 c o « N 09 V « N N 8 — 3 3cfi3«n ss a R RK R !S * <* - s 5 :I saps I s3 i fs ' ? ! s I M rt M j j i 5^ «A i T3 C a a: s 3 2 8 5 Mil 1111 s M * 8 ax a i- jl sl ? s SSI g iff o 8 SSI3 9« * ii 32Si ?5 E 11 I — r> -» §11 S R C• S 2 5§ M 3 8 asa s aaaa* » lis S g2 1 3aB 5 ii P SiSe 11 If: !§ ; R Appendix B Detailed Derivation of Pollaczek's Formula Fig. B.l shows the buried current filament in the earth. The properties of the air and the earth are given in the figure, x and y axes are established from origin O as shown in the figure, t axis is in the opposite direction of y. The location coordinates for the filament are (»/,!//) or (xf,tf) with tf = — yf > 0. region 1: air (M- ) fl O (x y ) r sf current filament • ,( x f f -I ,t ) f m \ Figure B.l: A current filament buried in the earth Applying the assumptions in Section 2.2 to this case the principal equations can be derived. Because of V x E = -juB (B.l) - V xB = J (B.2) V E =0 (B.3) 150 151 Appendix B. Detailed Derivation of Pollaczek's Formula the following equation is got VxVxE = V(V • E) — V E = - V E 2 = 2 -jwV x B = -jwpJ = ^ ^ E (B.4) P The current density has one direction only. Therefore, the two-dimensional principal equation descibing thefieldin the problem is V E = jvfiJ = \ E (B.5) P = -A- (B.6) 2 P where 2 p is the complex penetration depth. If E and E respectively represent Efieldsin the A E air and in the earth, (B.5) can now be splitted into two equations VE = 2 A VE 0 = 2 E \E. + jup.I6(x-x )6(t-t,) f Pi y>0 (B.7) t>0 (B.8) where Pe (B.9) "T~~ = jup, e I is the magnitude of thefilamentcurrent. The Delta function is associated with the filament. The boundary conditions are E =E =E a e ldE ldE a p, a ay a = y=0 0 ldE e e — = = p, ay dE — p, e a e = ~^ (B.10) = e dE e ~dT at dE y = 0; t = 0 dE a ~dV = = ^ ' e n = *= (B.ll) and/or y = oo (B.12) In order to solve (B.7) and (B.8) the integral transformation technique is applied. To x the following Fourier transformation pair is applied /(a) = / J—oo f{x)e~ *dx 3U (B.13) Appendix B. 152 Detailed Derivation of Pollaczek's Formula /(*) = ^ (B.14) J~J(a)e> da xa To y and t the following Fourier sine transformation pair is applied 71(A) = (B.15) [°° f(y)sm(Xy)dy Jo f(y) = - (B.16) rj {\)sm{y\)d\ t 7T JO For region 1, x and y are used as the coordinates. Define K = ja (B.17) E sin(Xy)dy (B.18) [°° E e- *dx a J—oo E = a / a Jo Applying the Fourier transformation to x as r°° dE 2 I J-oo dx e~ dx 3ax 2 = (-J«) = dEa dx - E e- (-ja) iax r / 2 / oo oo °° + f°° 2 a oo dE °° 8E -jax -oo E e~ dx }ax a E e- (-ja) dx jax = 2 a J-00 (B.19) -a E 2 a -00 d E~ 2 e~ dx = a 3ax By 2 (B.20) a dy 2 Then (B.7) becomes (B.21) Applying the Fourier sine transformation to y as roo J' (B.22) -a E s\n{Xy)dy = 2 a -a E 2 a f°° 0 1° d E2 . dy Jo svn\ Xy)dy 2 2 < dW dy . = -—— sin(Ay).oo dE dy a roo = - AjE cos(Ay)| - A y = XEQ — 2 o 0 XE 2 A E sm(\y)dy a (B.23) Appendix B. Detailed Derivation of Pollaczek's Formula 153 (B.21) is changed to T - ^ K (B.24) For region2, x and t are used as the coordinates. Define E r<x> / E e~ dx = e (B.25) 3ax e J — OO W e = rE~ sm(\t)dt Jo (B.26) e By similar process, the following equations can be derived -a T +^ dy 2 e {a + \ + \ )T 2 = ±--E +jupJe-i *'8{t-t ) pi (B.27) a 2 2 e f = \E '-jup Ie- ' m{\t ) (B.28) 0 =a +\ (B.29) jax e 0 e S f pi By assigning 2 2 Pi (B.28) becomes The inverse Fourier sine transformation is applied to (B.24) and (B.30) to get E and a (mi) = I* K According to the mathematical handbook f - ^ m Jo a* + M d X . * 2 a 7T S i n ( ; y > (B.34) y ! f <*A = JL ( « - ( « / - ' ) • _ 0 +A 40 ) 2 S (B.33) 0 y>0 2' r Jo > „ A) 2 ' v = i L ( -(t-t/)# _ -(*+t/)«) = J L ( -|t-t/ll«l _ -lW/ll«l) e e e e *>0; *,>*><> (B.35) 0>fj; *>*,>() (B.36) (B.37) Appendix B. Detailed Derivation of Pollaczek's Formula 154 Therefore, E a T. (B.38) = £oe" = E~ e-^ - juu.J -> *fJL |a|l/ ( -l«-*/ll«l _ -l-+-/ll«1) a 0 e e e (B.39) Using boundary condition (B.ll) E can be solved from the above two equations. Q Because dE a dy y=0 at = -\a\Eoe-M" y=0 = -\a\E (B.40) 0 t=o t=0 (B.41) the following equation is got - —\a\E^ = - — (-\9\EZ - ju, i.Ie-i 'e- 'W) am t t (B.42) Therefore, E is given by 0 ju>Ie- *e- tW jax EQ = jule- t t e ' \/ jax ±M + ±I'I -t a2+1 /p2 i H + i>/5 + i7ri (B.43) The final solution of 2? and E can be got by applying the inverse Fourier transfora e mation to (B.38) and (B.39) with E being replaced by (B.43). Also, t and tf can now 0 be changed back to — y and — y,, respectively. E and E are given by a 1 r°° E a = . —j 27T J — O O E^da v •°°iH + i \ / oo / / e e —/ ei( - t) da x a 2 + /Pe 1 2 x a ======= cos((x — z,)aWa (B.44) Appendix B. E Detailed Derivation of Pollaczek's Formula 155 1 r°° . = — / E e>* da a e e Jul r°° = -— / J 2TT = -*/N/« +i/rf-|«|t 2 e — = ^ - ^ da x x a 2|«9| 7-OO / , = cos((x — x f)a)da / 2TT y_oo , 2^/a + l/p 2 e>l -*A dct (B.45) m a 2 2? in the above equation can be further simplified into e E = e juu,I ( •i-fK (D/p ) - MD'IVe) 0 r°° 2e tW + /rt + I , (y+y e a2 ^ 1 cos((z - x,)a)da (B.46) where Ko(D/p ) = e / e>l —') da m a (B.47) •>-«> 2yJo? + I/pi roo -|-«-W/l\/ / j e Ko(D7Pe) = D a2+1 /pJ , Jl—'*da (B.49) = y/(x-x y + (y-y,)* = yf(x-xf¥ + (y + y y f f K is the zero order second kind modified Bessel funciton. 0 (B.48) (B.50) Appendix C List of Symbols A, A : magnetic vector potential [A] : A vector A : A function values on FEM boundary b [AB] : A B I A vector of Dirichlet boundary nodes element of [AB] : ak : distance to centre of PT cable of the A:th SC cable (&=A,B,C) a : unknown coefficient for trial function (p n O-mi • K : n polynomial coefficient in P (N ,() m p value of A at global node n A distribution in conductor k A (k) : A : A distribution in conductor k caused by conductor current J; : vector of A(ki) values in element Ei (LCI) [A^] [A ] = EI vector of A values at local nodes in element Ei element of [A ] Ei [Au] : A vector of unknown nodes B, B : magnetic field density [C] : C 0 i shunt capacitance matrix per unit length of transmission lines : coefficient in the formulas of SC coaxial cables C/,. : coefficient in the formulas of SC coaxial cables capacitance of thefcthinsulation in a SC cable 156 Appendix C. List of Symbols G : Ki Ct : coefficient in the formulas of SC coaxial cables diagonal element in [C] of SC coaxial cables C? : off diagonal element in [C] of SC coaxial cables C. = direct self capacitance of conductor i i0 G : mii Cik : direct mutual capacitance between conductors i and j element of [C] Cp : perturbation coefficient of the earth d : [D*] : differentiation shape function differentiation matrix in element Ei D : charge density D : distance between afieldpoint and a cable or diameter of strands D' : distance between afieldpoint and the image of a cable [D] : submatrix in [U] -f related to unknown nodes [Da] : submatrix in [U] + ju[T] related to boundary nodes [DL] = banded lower triangular matrix [Du] : banded upper triangular matrix common denominator in P (N , £) m p J - : width of the jth. division from a surface of the ith conductor tJ det[] : determinant of a matrix diag( ) : a diagonal matrix dl : integral element E, E : electrical field E : E value on the earth surface E : E field in the air E : E function values on FEM boundary 0 a b Appendix C. List of Symbols E B E c : vector of node values of E\, in a FEM mesh : Efieldin the earth with a deeply buried SC cable E. : E field in the earth E Efieldin the earth with a deeply buried current filament F : Ei : [Es] : E Si : element i in a FEM mesh source electricalfieldvector source electricalfieldin conductor t 70 : inverse Fourier transformation of /() [F] : integral coefficient matrix related to unknown nodes [F>] :update of [F] [FA] : a matrix modified from [F] [FB] = integral coefficient matrix related to Dirichlet nodes /« : division factor for the jth division Fk ' elements in matrix [F] FB,„ : elements in matrix [FB] Fi • F k in local node numbers in element Ei m E - mfe • m 1 Fs in local node numbers in element Ei lh [FH] = horizontal vector related to [F] IFk] ' update of [FB] TO • inverse Fourier sine transformation of /() [Fv] : a vector related to [EB] [G] : shunt conductance matrix per unit length of transmission lines [Gc] : diagonal conductivity matrix of conductors 9o(z,y) • Dirichlet boundary function h: burial depth of an underground cable or a current filament Appendix C. List of Symbols H : magnetic field h,2,fi3 : layout geometry parameters of tunnel installed SC cables i : index I : conductor current [7] : conductor current vector I , I eR ej I n : real and imaginary parts of the earth current : the 1st kind modified Bessel function of the nth order [IR], [II] : real part and imaginary part of [I] Io()j Ii() : first kind modified Bessel functions in 1st and 2nd order : conductor current in conductor i iRi hi '• real part and imaginary part of J; 5 iAi • internal return current of cylindrical conductor i iBi '• external return current of cylindrical conductor i 7f(r) : earth current within r of a deeply buried current filament [IEI], [IE2] '• equivalent current vectors due to boundary conditions [IE3], [IEA] '• equivalent current vectors due to boundary condition [EB] [1'EI], [I'E2] I ep I epc I ePF : updates of [I ] and [I ] E1 E2 : partial earth return current : partial earth return current of a deeply buried SC coaxial cable : partial earth return current of a deeply buried current filament j : complex number exclaimer or index J, J : current density «/(,-) : J distribution in conductor i J(ij) : J distribution in conductor i caused by conductor current Ij [ J(ij)i '• vector of J( ) values in element E\ tJ Appendix C. List of Symbols JR>JI '• 160 * l and imaginary part of J e a Js : source current density [Js] : source current density vector Jsi : source current density in conductor i [J ] : Jacobian transformation matrix in element Ei Ei k : index K : number of conductors K (), Ki() : second kind modified Bessel functions in 1st and 2nd order 0 Kf : factor in the GSW formula K n : the 2st kind modified Bessel function of the nth order [2/(u;)], [L] : series inductance matrix per unit length of transmission lines I : index Lc '. contour length of a sector-shaped core conductor [£c (<*>)] : series inductance matrix per unit length related to conductors [Lo] : series inductance matrix per unit length related to dielectrics Lij : elements of [L] M : number of elements in a FEM mesh m : index N : number of unknown nodes in a FEM mesh, number of space dimension n : index or the number of outer strands NB : number of Dirichlet boundary nodes in a FEM mesh NB : number of Dirichlet boundary nodes in the earth B : number of nodes in element Ei N P : degree of the auxiliary polynomials for simplex elements iVs : number of sampling points in an isoparametric element Appendix C. List of Symbols NT : total number of nodes in a FEM mesh NT = N + NB 0 : origin of x-y coordinates 0() : index sting for [QW] o(hi) : order of estimate errors p : complex penetration depth P : time-average power loss in a system or a field point Pm(Np,() '• auxiliary polynomials of degree N P p : earth complex penetration depth e Pij : time-average power loss related to conductor currents Ii and Ij Pi : complex penetration depth of conductor i Q : time-average reactive power in a system [q] : vector of surface charges on conductors qi : surface charge on conductor z qij : time-average reactive power related to conductor currents Ii and Ij Ag : surface charge on an element side : real integral matrices in simplex elements (k = 1,2,3) QW element of [Q™] : r : radius in polar coordinate system or the radius of outer strands r : FEM boundary radius b r : inner earth radius or the outer radius of a stranded conductor e r eq : equivalent radius from Borges da Silva's formula TA,TB '• inner and outer radii of an equivalent circular conductor *\Ai>rBj ' internal and external radii of cylindrical conductor i [i2(w)], [R] : series resistance matrix per unit length of transmission lines Rc ' internal resistance of stranded conductors from the GSW formula endix C. List of Symbols Rij : elements of [R] Re() : real part of the function s : separation distance of SC cables S : time-average complex power, or the eara of a triangle 5 ,5 : areas of subtriangles 2 3 Sc • cross-section area of a sector-shaped core conductor SR : solution region [Sc] ' conductor cross-section area matrix [S' ] • update of [Sc] c [Sc ] : a matrix modified from [Sc] A Sc : cross-section area of conductor k k [Sc ] : a vector related to [EB] v SEi '• region or area of element Ei t : time or t axis in the opposition direction of y axis tf : t coordinate of a buried filament [T] : coefficient matrix T mn : elements of [T] [T ] : imaginary integral matrix for isoparametric element Ei Ei : element of [T ] Ei [Ts] : imaginary integral matrix for simplex elements T Smn : elements of [T ] s [U] : coefficient matrix [U ] : [U] in element E Ei U mn t : elements of [U] Cfj, : U mn in local node numbers in element Ei Appendix C. List of Symbols u 163 : unit vector of z axis in Cartician coordinate system z u* : unit vector of 9 axis in Polar coordinate system [V] : conductor voltage vector with respect to the reference conductor Vi : conductor voltage of conductor i w : layout geometry parameter of tunnel installed SC cables WE ' electric energy found from the numericalfieldsolution WE : WE with conductors i and j energized P F F WE : electric energy from circuit analysis with conductors i and j energized C Wj : weighting factor of numerical integration for isoparametric elements WM : time-average magnetic energy stored in a system WMij '• time-average magnetic energy related to conductor currents Ii and Ij x,y : x and y axes perpendicular to the transmission line xi,X2,X3 : x coordinates of vertices in a simplex fiVi '• x x a n d y of a buried filament [x ] : vector of x coordinates of vertices for isoparametric element Ei Ei x n,y '• x and y coordinates of node n in a FEM mesh n xp,yp : x and y coordinates offieldpoint P yijjtejjfa : V coordinates of vertices in a simplex [y ] : vector of y coordinates of vertices for isoparametric element Ei Ei [Y(u>)], [Y] : shunt admittance matrix per unit length of transmission lines z : z axis parallel with the transmission line [Z(u;)], [Z] : series impedance matrix per unit length of transmission lines Z : earth return impedance e Zij : elements of [Z] Zki : submatrix in [Z] related to phases k and / (fc,/=A,B,C) 164 Appendix C. List of Symbols ZJH '• internal surface impedance of cylindrical conductor t Zsi ' external surface impedance of cylindrical conductor i Z\n '• transfer impedance of cylindrical conductor i Zf. : diagonal element (k, k) in [Z] of a SC coaxial cable Z%? : off diagonal element (i, k) and (k,j) < k) in [Z] of a SC coaxial cable ZsQi : equivalent impedance related to cylindrical conductors i nad i + 1 Zj) : impedance related to the dielectrics between cylindrical conductors i and j i [ly] : a vector filled with 1 and 0 [ ] : transposition of the matrix £ [ ]* : conjugate of the function | | : absolute value of the function a : variable in Fourier transformation mimm '• shape function for simplex elements : 2 3 [/?] : shape function vector for isoparametric elements Pi : shape function for isoparametric elements S : real penetration depth in conductors 6i : S in conductor i 8 : 8 in the earth e £() : Deric function 8 : real penetration depth in conductors e : permittivity €Ei '• permittivity in element Ei Cfc : permittivity of the k insulation c : relative permittivity r £ : simplex coordinate Appendix C. List of Symbols d • simplex coordinate i 9 : a variable or span angle of sector-shaped regions #: include angle of vertex k in simplex element Ei A : wave length or variable in Fourier sine transformation V'o : function to enforce the Dirichlet boundary condition fl : permeability Ha • permeability in the air fi : e permeability in the earth /*o : permeability in the vacuum fi : r relative permeability A*fc : permeability in conductor k /** : /> : : permeability in element Ei volume charge density or conductor resistivity earth resistivity <7 : conductivity <Tk •conductivity in conductor k O'Ei ' conductivity in element Ei : local coordinates for isoparametric elements : sampling point local coordinates in isoparametric elements : local coordinates for isoparametric elements electrical scalar potential or position angle electrical scalar potential vector value of <f> at global node n 4* : local node value of d> in element Ei [d>] in element Ei Appendix C. List of Symbols u> : radius frequency <Pn : VBi : trial function for global node n trial function for Dirichlet node t : trial function for local node n in element Ei r : boundary surrounding the solution region r : Dirichlet boundary 0 Ti : rc> : homogeneous Neumann boundary periphery of the cross-section area of conductor
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Calculation of frequency-dependent parameters of underground...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Calculation of frequency-dependent parameters of underground power cables with finite element method Yin, Yanan 1990
pdf
Page Metadata
Item Metadata
Title | Calculation of frequency-dependent parameters of underground power cables with finite element method |
Creator |
Yin, Yanan |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | In this thesis, the finite element method (FEM) is applied to the calculation of frequency-dependent series impedances and shunt capacitances of underground power cables. The principal equations describing the quasi-magnetic fields and static electric fields are solved with FEM based on the Galerkin technique. The Js method and the loss-energy method are derived to calculate the impedances of a multiconductor system from its field solution, and the energy method and the surface charge method are derived to calculate the capacitances. With a single-core (SC) coaxial cable, the suitability of quadratic isoparametric elements and high-order simplex elements are studied, and a suitable division scheme is suggested for the auto-mesh program. The conventional FEM with a field truncation boundary is applied to the impedance calculation of buried SC cables. Suitable locations for the field truncation boundary and division schemes in the earth are studied. The results show that rb ≥ 12[symbol omitted] is required to obtain accurate impedances of shallowly buried cables with the conventional FEM. This requires a large solution region in the earth at low frequencies. A new technique based on the perturbation concept is proposed to reduce the solution region in the earth. Comparisons between the results from the conventional FEM and from the proposed technique with a significantly reduced solution region in the earth show good agreement. In the case studies, the FEM is applied to the parameter calculation of multiphase SC cables, PT cables, sector-shaped cables, and stranded conductors. The numerical results are compared with those from analytical formulas. |
Subject |
Underground electric lines Finite element method |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-08 |
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.0100550 |
URI | http://hdl.handle.net/2429/31119 |
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 |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1990_A1 Y56.pdf [ 11.32MB ]
- Metadata
- JSON: 831-1.0100550.json
- JSON-LD: 831-1.0100550-ld.json
- RDF/XML (Pretty): 831-1.0100550-rdf.xml
- RDF/JSON: 831-1.0100550-rdf.json
- Turtle: 831-1.0100550-turtle.txt
- N-Triples: 831-1.0100550-rdf-ntriples.txt
- Original Record: 831-1.0100550-source.json
- Full Text
- 831-1.0100550-fulltext.txt
- Citation
- 831-1.0100550.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-0100550/manifest