CALCULATION OF FREQUENCY-DEPENDENT PARAMETERS OF UNDERGROUND POWER CABLES WITH FINITE E L E M E N T METHOD By YANAN YIN B. Eng., Shanghai Jiaotong University, 1982 M . Sc., Tsinghua University, 1984 A THESIS 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 OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F 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 agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia Vancouver, Canada Date StifMtvtU^ 26j m o DE-6 (2/88) Abstract 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 capac-itances. 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^ , > 12Se 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. 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 [Z] Calculations from the Field Solutions 20 2.5.1 The Js method 20 2.5.2 The loss-energy method 21 2.6 Impedances of Single-Core Coaxial Cables 24 2.7 Summary 26 3 Element Types and Shape Functions 28 iii 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 Isoparametric Elements 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 42 3.5 General Procedures for [Z] calculations with FEM 42 3.6 [Z] Calculation of an SC Coaxial Cable with FEM 45 3.6.1 Optimum divisions for SC coaxial cables 47 3.6.2 Computation efficiency: CPU time, storage, and pivoting . . . . 55 3.7 Summary 58 4 Earth Region Reduction Technique for [Z] Calculation with F E M 60 4.1 Introduction 60 4.2 Analytical Formulas for [Z] of Buried SC Coaxial Cables 62 4.3 [Z] Calculation of Deeply Buried SC Coaxial Cables by Conventional FEM with a Field Truncation Boundary 67 4.4 Earth Reduction Technique 74 4.5 [Z] Calculations of Shallowly Buried SC Coaxial Cables with FEM . . . 82 4.5.1 Determination of the solution region for the conventional FEM in [Z] calculations of shallowly buried cables 83 4.5.2 Application of the proposed technique to [Z] calculations of shal-lowly buried SC cables 86 iv 4.5.3 Comparisons between analytical results and FEM results for shal-lowly buried SC cables 91 4.5.4 [Z] calculations for a cable layout of arbitrary structure 93 4.6 Summary 95 5 Admittance Calculation with Finite Element Method 97 5.1 Introduction 97 5.2 Principal Equation and FEM solution 98 5.3 [C] Calculation from the Field Solutions 100 5.3.1 The energy method 101 5.3.2 The surface charge method 103 5.4 [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 5.5 Summary 112 6 Case Studies in [Z] and [C] Calculations 114 6.1 Introduction 114 6.2 [Z] Calculations of Buried or Tunnel Installed Multiphase Cable Systems 115 6.3 [Z] and [C] Calculations of PT Cables 119 6.3.1 The [Z] calculation of PT cables 119 6.3.2 The [C] calculation of PT cables 125 6.4 [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 137 References 141 A Integral matrices [Q^] and [Ts] of simplex elements 146 B Detailed Derivation of Pollaczek's Formula 150 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 10n to 10n + 1 70 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 . . . . 79 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 4.8 Iep found with numerical integrations for the proposed technique . . . . 87 4.9 Storage and other parameters for the proposed technique (rfc=5m) . . . . 89 4.10 CPU time requirements for the proposed technique (rt=5m) 89 4.11 Storage and other parameters for the conventional FEM (rj,=12£e) . . . . 89 4.12 CPU time requirements for the conventional FEM (rj,=12£e) 90 4.13 Maximum differences in [Z] with the proposed technique at different pe . 90 4.14 Maximum differences in [Z] with the proposed technique at different r e . 90 4.15 Maximum differences in [Z] with Pollaczek's formula at different pe . . . 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.19 [Z] of the tunnel installed cable from three approaches 94 5.1 i m ando m i inP m (A r p , 0 = i E f c 0 o m i C < 106 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 I l l 6.1 [ ^ A A K [ ^ C C D A N < ^ [ ^ B B I °^ D u ried three-phase cable system 117 6.2 [ • Z A B K ^ B C D a n < ^ [^ACl °^ buried three-phase cable system 117 6.3 [ ^ A A K I ^ C C D a n d [^ BB] °^ * n e tunnel installed three-phase cable system 118 6.4 [ • ^ A B K ^ B C ] ) [^ACl °^ * u n n e ^ installed three-phase cable system 118 6.5 [ ^ A A I A N ( ^ [ ^ B B K ^ C C } ) °^ c a ^ e with a triangle arrangement . 123 6.6 [ • Z A B K ^ A C D a n < ^ [^ BCl °f * n e cable with a triangle arrangement . 123 6.7 [ ^ A A K I ^ C C I ) [ ^ B B I °^ PT cable with a cradle arrangement . . 126 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) 127 viii 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 146 A.2 [Q^\ and [Ts] of the 2nd order simplex element 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 [Ts] of the 5th order simplex element 148 A.6 and [Ts] of the 6th order simplex element 149 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 (JVp=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 Program flow chart for [Z] calculations with FEM 46 3.9 Geometry of a SC coaxial cable and its FEM solution region 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 A shallowly buried SC coaxial cable x 63 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 filament 77 4.6 Perturbation coefficient Cp as a function of |r e/p e | 78 4.7 FEM mesh at 6kHz for rb=2Se 83 4.8 Earth return current and maximum differences in [Z] for different r& . . . 85 4.9 E field at 6kHz given by the analytical solution 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 under-take 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 frequency-dependent 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 core air earth SC coaxial cables (b) Tunnel-installed SC coaxial cables (a) Different types of three-phase 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 electromagnetic fields inside 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 pa-rameters are then derived from the numerical field solution. 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 prox-imity effects, it is more difficult to solve quasi-static magnetic fields than 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 was first applied to the series impedance calculation of multicon-ductor 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 very fine subdivi-sion has to be used to achieve accurate results if strong skin and proximity effects are present. With the finite element method, the domain of a closed boundary problem is divided into small elements such that the original unknown field distribution in each element can be approximated by certain functions expressed in terms of unknown field variables at element vertices. The field variable 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 to find the field distribution 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. The final linear 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 coeffi-cients 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 prob-lems [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. The final coefficient 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 develop-ment 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] calcu-lations, 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 con-ductive 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 under-ground power cables, i.e. Laplace equations, are solved with the FEM, and the formula-tions 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 admit-tance 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 Chapter 2. Impedance Calculation with Finite Element Method 9 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 cP, permeability /zP, and conductivity a. 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 between field quantities, 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\, V2, ..., VR are the conductor voltages with respect to the reference conductor K+l, and / j , 72, . . . , IK 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_- -Jg^ , conductor k . {*+«'* Figure 2.1: A (iif+l)-conductor system [/] [LD)dz [R((a)]dz {1^(0!)] dz [V] O [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 M dz = ([R(u,)}+ML(u)])[I] = [Z(u)][I\ = ({G]+MC])[V) = [Y(u)}[V} (2.1) (2.2) in which [V] =[VUV2,...,VK? (2.3) Chapter 2. Impedance Calculation with Finite Element Method 11 [I] =[IUI2,...,IK]T (2.4) [Z(w)] = [R(v)] + ML(u)] = [J2(w)] + M[Lc(u>)] + [LD]) (2.5) [*»] =[G]+MC] (2-6) [Z(u;)] and [Y"(w)] are the series impedance matrix per unit length and shunt admit-tance 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 magnetic field and a static electric field of 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 magnetic fields have 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 assump-tion 4. As a result, the original electromagnetic field becomes 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, Wede-pohl and Efthymiadis rigorously analyzed the overhead line case over the full frequency Chapter 2. Impedance Calculation with Finite Element Method 12 spectrum[18][19]. They found that the assumption was valid until the height of the con-ductor 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 be-tween field tests 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 = V x A (2.10) and inserting (2.10) into (2.8) and employing unity V x V x A = V(V • A) — V 2 A , the following equation is obtained - V ( V • A) - - V 2 A = J (2.11) ft ft Assuming V A = 0 (2.12) (2.11) gives - i v 2 A = J (2.13) As the current density has only a longitudinal component, J , E, and A can be written respectively as Ju z , Euz, and AVLZ. U z is the unit vector along the z axis. By inserting Chapter 2. Impedance Calculation with Finite Element Method 13 (2.10) into (2.7), the following equation can be derived (E + jwA)uz = -V<f> (2.14) 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 s u z = (rEsuz = —<rVd> (2.15) (2.14) can now be written as J = -jvcrA + Js (2.16) Combining (2.13) with (2.16), the following linear two-dimensional diffusion equation can be derived - V 2 A - jwvA + Js = 0 (2.17) The integration of (2.16) over the cross section of a conductor gives 1=1 Jds = -ju f a Ads + ScJs (2.18) JSc 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. Chapter 2. Impedance Calculation with Finite Element Method 14 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] cal-culation can be summarized as -V2A-ju(rA + JS = 0 inSR (2.19) A* -Jul <7Ads + SCkJsk = h {k = l,2,...,K) (2.20) JsCk with boundary conditions = 0 (2.22) ri ^lr0 = 9o(x,y) (2.21) dA dn SR is the solution region surrounded by boundary r = ro + IV Sck and Jsk are the cross-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 the field distribution A satisfying the differential equation (2.19) by a finite set of base functions ipn (n = 1,2,..., N) as N A = ip0 + Y, °n<pn (2.23) where ax, a 2 , . . . , ajv & i e unknown coefficients. ipQ is such that Mv0=% (2.24) V'lj ¥>2> • • -i VJV form a subset of the complete base functions of A. ipn satisfies the boundary condition Vn|ro=0 (n = l,2,...,JV) (2.25) Chapter 2. Impedance Calculation with Finite Element Method 15 Putting the approximate solution (2.23) into (2.19), a residual is introduced as R(A) = -V2A - juaA + Js in SR (2.26) /* The Galerkin technique forces this residual to satisfy the following integral equation / R(A)ipnds = 0 (n = 1,2,...,iV) (2.27) JsR 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 + £ an<pn) \ <pmds - f jW(^0 + £ anipn)<pmds + JsR l*\ £i J JsR + / Js<Pmds = Q (m = l,2,...,AT) (2.28) JsR Applying Green's formula V • (vVu) = Vu • Vu + vV2u, the first term in the above equation becomes r 1 N r 1 N - / -Vy?m • V(V>o + £ an<pn)ds + - V • (y>mV( 0^ + £ anVn))ds JSRfi n = i JSRH N = 1 = - / - V ^ - V ^ o + E a n ^ + Z I ^ m ( ^ + £ a J ^ r JSR^ n=i Jr°+rl I1 9 a n=l 9 n (m = 1,2,..., AT) (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 <pn Chapter 2. Impedance Calculation with Finite Element Method 16 is established, the N unknown coefficients are found by solving the following equations 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 equa-tions (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 An (n = 1,2,..., N), with N being the number of nodes within the mesh. The coordinates of An are represented by (xn, yn). The shape function y>n (n = 1,2,..., N) in the [Z] calculation satisfies the following conditions. (m = l,2,...,iV) (2.30) Chapter 2. Impedance Calculation with Finite Element Method 17 1. ipn is continuous inside SR. Therefore, it is continuous across the interelement 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 0 with prescribed values. A similar expression can be written for ipo as NB ^0 = YLABi lPBi (2.31) t=l in which AB% represents the known node value on the boundary To, <fiBi represents the corresponding shape function which has the same characteristics as <pn discussed above, 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 A = Y£AN<PN (2.32) n=l and (2.30) and (2.20) become r 1 N t r N t r - / -V<^m • V J* AnVnds - I j<jJ<T(pm J2 Anipnds + / Js<pmds JsRfi ^ JsR ^ JsR = - V) AN I (-V<pm • V<pn + ju}cripm(p^\ ds+ I JsVmds = 0 (m = l,2,...,iV) (2.33) Chapter 2. Impedance Calculation with Finite Element Method 18 - J ^ M I <7<Pnds) + SCkJsk = h (k = 1,2,..., A") (2.34) n=l "/SCfc Equations (2.33) and (2.34) are combined together to form the following final algebraic equations in matrix form [U]+ju[T] \ -MGo] [ [F]T, [FBY ] [So] [0] (2.35) in which m = 1,2,....W n = 1,2,...,NT l = N + l,N + 2,...,NT Jb = l,2,...,JT (2.36) Umn = ISR ^V< m^ • V<pnds Tmn — fsR VWmVnds F M K = ISCk fmds FB* = IsCk fids [Gc] = diagfa,a2,... ,crK] [Sc] = diaglScnSct^.-yScK] [A] = [AUA2,...,ANT)T [Js] = Psi > Js2) • • • j ^5K]r 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]T, [FB]T ]• In practice, (2.35) is assembled element by element. That means all the integrals in (2.36) become the summations of integrals over the elements as [ [ (2 3 7 ) JSR JSEi 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 SEM , and Sc* £ SR. Inside the ith Chapter 2. Impedance Calculation with Finite Element Method 19 element (2.32) becomes NEI A=E A*<p« in SEi (2.38) n=l in which NEi is the number of nodes in the ith. element, AEi is the node value of A, and ipEi is the shape function defined locally in the element. Umn and T m n become UEi = • V(pEids rpEi _ mn JSBi ds (2.39) J?Ei — JsEi for conductor elements *%> = f <fiffds for conductor elements where m = 1,2,...,NEi excluding boundary nodes, n = 1,2,..., NEi, I = 1,2,..., A ;^. excluding unknown nodes, and i = 1,2,..., M. Boundary nodes are those on the bound-ary T 0 with prescribed node values whose global numbers are bigger than N, and unknown 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 fiEi and o~Ei are the permeability and conductivity in the 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 magnetic field represented by the discrete node values of A can be found by solving (2.35). From the field solutions, the series impedance can be calculated. This is the topic of the next section. Chapter 2. Impedance Calculation with Finite Element Method 20 2.5 [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 the field solution 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 Js method As mentioned in Section 2.3, the physical meaning of Es in (2.15) from the field analysis 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 [Es] = [ESi,ES2, • • •, ESK]T- Therefore, [Js] = [GC][E8] = -[Gc]^- = [Gc][Z][I] (2.41) If [/] is assumed as Ii = 0 and Ij ^ 0 = 1,2,..., K; i ^ j), the jth column of [Z] can be derived from the above equation as Zii = Ti- (« = 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 time-average reactive power in the equivalent network. Assuming [I] = [IR] + j [2j] and using the fact that P = PT and Q = QT, the following equations can be derived P = [IR]T[R}[IR] + [II]T[R][II} = TiTIPV (2-44) t=l j=l K K Q = o,([//[i][/fl] + [/;f[I][//]) = E E ^ (2-45) t=i i=i in which Pij = RiiWRi + hh) (2-46) Q I I = UL^IRJR, + //.//.) (2.47) The time-average magnetic energy WM stored in the system is related to Q by WM = ±Q = hlR}T[L][IR] + [//]T[I][//]) = E E wMii (2.48) in which wiiii^LiiWRi+hhi) (2-49) Consequently, if 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* p = Y.L ~irds (2-5°) fc=l JSck <* WM = I I BH*ds = lT/Re( [ AJ*ds\ (2.51) 2JsR 2 £ i \Jsck J 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(K) is the current density in the kth. conductor, it can be written as K J(K)=T,J(KI) (A: = 1,2,..., it) (2.52) t=i where is the current density in the kth. conductor caused by the current /,• in the ith. conductor. (2.50) now becomes k=l Jsck a k=\ Jsck a i=i ,=i K . -i K K K K K , 7.,.. = E / s ; E E W = E E E / S ™ * k=l JSCk a <=1 j=\ i=l 7=1 k=l JSC„ a 3-K K = EE ^J (2-53) »=ii=i P I J — J* k=lJbck = E / ™ ^ (2.54) k=lJsck a Similarly, = £ £ o £ R e U = E E ^ (2.55) .'=1 j=l Z A=l V-76^ / «=1 j=l ^ = ^ E R e ( ^ c V-) J (V 5 ) (2-56) where = J2iLi ^.(Ai) is the A in the Mh conductor caused by current /,-. Relating (2.54) and (2.56) with (2.46) and (2.49), respectively, gives and L{}- as **i = {EL Jjt^ds)l{IJUIRi+IIiII}) {i,j = l,2,...,K) (2.57) Kk=lJSck a J La = ^«)^y)*)}/(J^J*i+V/i) (i,i = l,2,...,A')(2.58) Chapter 2. Impedance Calculation with Finite Element Method 24 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 ana-lytically 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 the fcth cylindrical 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 electric field Euz satisfy the following equations in cylindrical coordinates E=L4±m ( 2 . 6 0 ) ro~ dr The solution for E inside the Arth conductor is E{r) = — {CIkh(r/Pk) - CKkK0(r/pk)) rAk < r < rBk (2.61) o~Pk where CIk = 7T7^{—Ki(rWP*) + —KxCrxJp*)) ZTrL>£)k \rAk rBk ) CKH = « ^ f^Ii(rB*/w) + ^ IiCrxik/p*)) (2-62) ZnCDk \rAk rBk ) CDk = h(rBk/Pk)Ki(rAk/Pk)-IiirAk/P^Kxirek/pk) Pk = \ -r^ V JUCTkHk IAk and IBk a r e internal and external return currents of the fcth conductor, respectively. pk is a complex penetration depth in the kth conductor. In, and K n (n = 0,1) are the modified Bessel functions of the nth order and respectively of the first and second kinds. E on the internal and external surfaces of the fcth conductor can be derived as E(rAk) = ZAkIAk + ZMkIBk (2.63) E(rBk) = ZMkIAk + ZBkIBk (2.64) where zAk = 2TrrAklkpkCDk ( I °( 7 '^ f e /^) K l ( r f l fe/P f c ) + Ko(rAk/Pk)li(rBk/Pk)) (2.65) zBk = 2irrBkakpkCDk (^( r B f e/P f c)K l( r^/P f e) + Ko(rBk/Pk)h(rAk/pk)) (2.66) 1 , Chapter 2. Impedance Calculation with Finite Element Method 26 ZAICI Zsk* a n < i %Mk a r e defined by Schelkunoff as internal surface impedance, external surface impedance, and transfer impedance of the fcth cylindrical conductor, respectively. Suppose that the system in Fig. 2.4 has a superconductive surface at TAK+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] where ZEQ K t [Z] Zd yod yod 1 L7. "3 yod yd yod ^2 L1 L3 \ Zod yod 3 ^ 3 zt = zf JEQk = yod yod yod \ LK LK LK K K Yl ZEQK - 2 E ZMk k-i k=i+l K K ZzQ-k ~ ZMX - 2 ]^ Z\fk %Bk + %Dk + %Ak+l z$ z$ z$ Zjc J = ZBK + Z2 (i = l,2,...,K) (i = 2,3,..., K) (k = l ,2 , . . . , t f - l ) (i = l,2,...,K) (2.68) (2.69) (2.70) (2.71) (2.72) (2.73) 2TT \ rBk , ZAk+ii ZBHI and Z\ik a r e 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 tech-nique. The Js method and the loss-energy method are derived for the impedance calcu-lation 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>E' (n = 1,2,..., A T E J inside an element is still unknown. The shape function <f>Ei is related to element shapes and interpolation functions in the elements. In 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 fre-quencies. 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 curve-sided 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] calcu-lations. 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 the field solutions are compared. 3.2 High-Order Simplex Elements The approximate field solution given by (2.32) in the preceding chapter can be viewed as a two-dimensional interpolation formula. Once the field values at discrete nodes are known, the field value 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 in finite element analysis, because they can be differentiated and integrated easily. The corresponding shape functions of the simplex elements and the isoparamet-ric 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 polyno-mials 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 Chapter 3. Element Types and Shape Functions 30 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 the figure. 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 coordi-nate system. Only two of the three simplex coordinates are independent because Ci+e2+c3 = 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). Chapter 3. Element Types and Shape Functions 31 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 ( 322/3 - x3y2 2/2 - 2 / 3 « 3 - « 2 \ 1 = 25 332/1 - « l J / 3 2/3 - 2 / i xx - x3 X ) ^ Z l l / 2 - Z2J/1 2/i - 2 / 2 x2 - 3 i / \y ) In a simplex the shape functions can be expressed in terms of simplex coordinates as [34] (Np,Ci)Pm2(NP,C2)Pm3(NpiC3) m1+m2+m3 = Np (3.4) in which Pmi(Np,Ci), Pm2(Np,C2), and Pm3(Np,£3) are auxiliary polynomials given by Pm(Np,C) = A l l W - f c ) m = l,2,...,AT p (3.5) m ! fc=0 Po(ATp,C) = 1 (3.6) Np is the order of the polynomials. Integer indices n»i, m2, and m 3 have values from 0 to Np and are related to the locations of nodes inside the triangle. They will be 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 Np 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 Chapter 3. Element Types and Shape Functions 32 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 (iVp=4) 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^, ipEi = a(n_!)10,..., and ip^ = a0on' For the nth order 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%n in (2.39) can be derived as[34] = 7 - E (3-7) Chapter 3. Element Types and Shape Functions 33 NEl-l Figure 3.3: A commonly used node numbering scheme (3.8) in which V m n Js* \dCk+l dCk-x) \d(k+i dCk-iJ 2SEi JSA, &E: (k = 1,2,3) (3.9) (3.10) 0E< 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 d3 = dxdy = dx 8x By By S<1 B<2 d(id£2 = 2SEidCid£2 (3.11) 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 r1 r1-^ . . . . . . + * + 2)! (3.12) Chapter 3. Element Types and Shape Functions 34 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 Order Index string 0() l 3 1 2 2 6 3 5 1 2 4 3 10 6 9 3 5 8 1 2 4 7 4 15 10 14 6 9 13 3 5 8 12 1 2 4 7 11 5 21 15 20 10 14 19 6 9 13 18 3 5 8 12 17 1 2 4 7 11 16 6 28 21 27 15 20 26 10 14 19 25 6 9 13 18 24 3 5 8 12 17 23 1 2 4 7 11 16 22 F^k and F§*k can be calculated from [Ts] as * S = SEiJ2Tsmn (3.13) n=l = (3-14) n=l where m = 1,2,... ,NET, excluding boundary nodes, and I = 1,2,... ,NE{, excluding 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 in-terpolation 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) Chapter 3. Element Types and Shape Functions 37 © f ^ ^ D , E. (a) global coordinates (-1.D © © © (-1.-D © (-1.0) (0.1) (1.1) ' © (1.0) © © 1) © (0.-1) (1,-1) (b) local coordinates Figure 3.5: Quadratic quadrilateral isoparametric element x = \j3][*Ei) y = \P}[yEi) (3.16) (3.17) where [AE<] [**] — F l > • c2 > • ..,*?]* far*] = [yfSyfv (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/) f32 = \{l + v)(l-v)(-l+v-v) P3 = i(l + v)(l + u)(-l + v + u) P4 = \(l-v)(l + v)(-l-v + v) & = 1(1 - v2)(l - *) fh = i(l-v'){l + v) /?8= 1(1-^(1-^) (3.19) Chapter 3. Element Types and Shape Functions 38 From (3.16)-(3.18), the following equations can be derived a a dx = [J?]"1 dv a a L ay J . dv . ds = dxdy = det[JE']dvdv where [JEi] is the Jacobian matrix of the transformation given by (3.20) (3.21) -ax fit " a dv av dv dx a . dv dv . W ([**] [yEi]) (3.22) The integrals U$n and T$n in (2.39) can now be evaluated in the local coordinates. The matrix form will be used. [UEi] and [TEi] are the matrices corresponding to and TE'n, respectively. These matrices are 8x8 square matrices. Some of their elements, however, may not be used in the final matrix assembly because subscript m in U^n and TE^ only refers to the unknown nodes in the element. [UEi] and [TEi] are given by ^ = V~L V[/3] • V[/3]r^ = / UEi JSBi f*Ei JSBi = — f1 f1 [DEi)T[DEi}det[JEi]dvdv fiEi J-i J-i [TEi] = <TEJ^JP]T[P}det[JE<}dvdv a T ( a dx dx a \ a I 8y J I L Sy J [/3] ds in which [DE<] = a a dx Ifl = [Jf]"1 dv a a I ay J . 8v . [/?] (3.23) (3.24) (3.25) As the integrals in (3.23) and (3.24) are related to the element node coordinates [xEi] and [yEi], [UEi] and [TEi] cannot be evaluated once and for all. Also it would be too complicated to integrate (3.23) and (3.24) analytically. Instead, numerical integration Chapter 3. Element Types and Shape Functions 39 formulas are used. With Gaussian quadrature formulas, [UEi] and [TEi] become 1 Ns Ns [UEi] = — E E ^ [ ^ * ( * i . ' * ) ] r [ ^ * ( « i ^ ) ] M ^ , ( « i ^ ) ] (3-26) VEi j=l k=l Ns Ns [TEi] = ^ E E ^ ^ ^ > ^ ) ] r ^ i , ^ ) ] ^ [ J f ( ^ , ^ ) ] (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 quadra-ture 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), FEih and FE* are calculated from [TEi]. Table 3.2: Locations of sampling points and weighting factors for Gaussian quadrature Ns Vj or Vj Wj i 0 2 2 1 3 0 8 9 ±0.7745966692414834 5 9 4 ±0.3399810435848563 ±0.8611363115940526 0.6521451548625461 0.3478548451374539 5 0 ±0.5384693101056831 ±0.9061798459386640 0.5688888888888889 0.4786286704993665 0.2369268850561891 3.3.2 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 (a) global coordinates (0.0) {05,0) (1,0) V (b) local 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 = a2oo = Ci(2& - 1) = v(2v - 1) 02 = "020 = C2(2C2 - 1) = v(2v - 1) 03 = «oo2 = Ca(2C3 - 1) = (1 - v - -2v- 2v) (3.28) 04 = "no = 4CiC2 = 4w 0s = "on = 4C2C3 = 4i/(l - v - v) 06 = aioi = 4CiC3 = 4v(l - v - u) The numerical integration formula is applied again to calculate [UEi] and [TEi] \UEi\ = ^ E ^ i l ^ ( « i ^ i ) ] T P * ( « i » ^ ) ] ^ [ ^ ( « i . ^ M (3-29) [TEi] = ^ E ^ ^ i ^ i ) ] r ^ i ^ i ) ] ^ [ ^ K ^ i ) ] (3-30) 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 Chapter 3. Element Types and Shape Functions 41 polynomial of order k — 1. Compared with simplex elements, fewer elements are needed if the above isoparamet-ric 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 isoparamet-ric 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. FEik and FE*h are calculated from [TEi]. Table 3.3: Locations of sampling points and weighting factors for quadratic triangular element Ns Wi Error i (3 ' 3) 1 2 0(^2) (0,0.5) 1 6 3 (0.5,0) (0.5,0.5) 1 6 1 6 o(h3) 4 V3> 3/ (11 J.) M5> 15 / (2. 11) M5> 15/ (± ±) V15> 15/ 27 "96 25 96 25 96 25 96 o(h4) V3> 3 / 0.1125 (0.47014206,0.05961587) 0.066197075 (0.47014206,0.47014206) 0.066197075 o(h6) 7 (0.05961587,0.47014206) (0.10128651,0.10128651) (0.79742699,0.10128651) (0.10128651,0.79742699) 0.066197075 0.06296959 0.06296959 0.06296959 Chapter 3. Element Types and Shape Functions 42 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 J j ^ J f c ^ = £ ^ • nT,][J*J] (3-31) in which I refers to elements in conductor k, [J^] and [Jr y^)] are node value vectors of current density distributions «/(«) and Jfa in element Ei, respectively. is the conjugate of J(kj)- Similarly, IOM,-,- in (2.56) is given by wMij = E Re / Awfads) = E Re fe SE,[A^]T[TS)[J^]\ (3.32) k=l \ l JSBI / 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, and WM^ are respectively given by = E E A f t ] r ^ ' ] f t l (3-33) fe=l l aEi wMij = E R e f E ^ - [ 4 o ] r ^ ' ] [ « 7 S ) ] ) (3-34) k=i \ i Et ) where [TEi] is given either by (3.27) for quadrilateral elements or by (3.30) for triangular 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). Chapter 3. Element Types and Shape Functions 43 As mentioned in Section 2.4.2 that boundary nodes are numbered after unknown nodes, vector [A] can be partitioned as [A] = [Av] [As] where [Av] = [Ax, A2,..., AN]T and [AB] = [AN+1,AN+2,, [U] + jw[T] is partitioned as [U]+MT} = [D] [DB] (3.35) ,, ANT]T. Correspondingly, (3.36) where [D] is an N x N sparse complex symmetric matrix, [DB] is an N x NB complex matrix. With (3.35) and (3.36), (2.35) becomes in which (3.37) ' [D] ~[F}' ' [Av] ' [hi] [Sc] . . m + [ / « ] . [FB] = -MGC][F]T [IEI] = -[DB][AB] [ISI] = MGC][FB]T[AB] (3.38) [IEI] and [7^ 2] are equivalent current vectors due to Dirichlet boundary conditions. In most cases of [Z] calculations, the boundary value vector [AB] is a zero vector. 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] = [DL][Dv] (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 (3.40) in which [DL][F'\ = —[F] [DL}[I'EI] = [IEI] 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] ' . [FH] [So]m . M . . [I] + [/«] . (3.41) '[Du] [F>] ' ' [Au] ' WEI] . PI [S'c]. . [Js) . . [I] + [Ik*] . (3.42) in which [S'c] [I'E2\ [F'H][Du] [Sc]-[F'H][F'] [lE2]-[F'H][I'E1] [FH] (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 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. [S'c][Js] [Du][Au] [I] + [iiJ (3.44) Chapter 3. Element Types and Shape Functions 46 (Start) / Read data / Process for simplex elements I Renumber nodes 1 Form [U], [T], IF], [FB], and [Gc] Assign / Form [[D][DB]] = [U]+j(o[T] andget [[ D ] Aufl _ I" [ /£1] "1 IFHI [ 5 c ] J | [ / , ] J " [ t / M / £ 2 ] J Factorize [D] = [D u][DL ] I Assign [/ ] Solve [/ s] and U J Calculate [Z] /"Plot fields and print [Z ] / Yes (Stop) Figure 3.8: Program flow chart for [Z] calculations with FEM Chapter 3. Element Types and Shape Functions 47 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. The first step 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. The fineness of the mesh will affect not only the accuracy but also the amount of computation. In general, a fine mesh 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 the fineness of 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 conduc-tivity, 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 Chapter 3. Element Types and Shape Functions 49 written as da = fdiSi (3.46) 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 Si. /a. siml 0.3 0.5 0.9 1.8 iso or sim2 1.15 2.8 — — sim3 2.25 — — — sim4 3.1 — — — sim5 4.1 — — — sim6 5 — — — "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) division radii (mm) 6 0 12 18 22 24 60 0 12 18 22 24 600 0 8.87 12 18 22 24 6000 0 8.601 11.01 12 18 20 22 24 60000 0 10.925 11.687 12 18 19.079 20.921 22 24 600000 0 11.66 11.901 12 18 18.341 19.171 20.829 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) L (jiH/km) /(Hz) element Rn R\2 R22 i n L22 ana 0.0388114 0.216122x10" B 0.414466 188.610 36.1306 29.4773 iso 0.0388114 0.215176x10" • 6 0.414466 188.590 36.1314 29.4758 siml 0.0388839 0.215349x10" • 6 0.415225 187.280 36.1247 29.2811 sim2 0.0388823 0.214880x10" 6 0.415225 188.592 36.1261 29.4689 6 sim3 0.0388823 0.215736x10" 6 0.415225 188.607 36.1277 29.4742 sim4 0.0388823 0.215739x10" 6 0.415225 188.607 36.1282 29.4747 sim5 0.0388823 0.215740x10--6 0.415225 188.607 36.1283 29.4748 sim6 0.0388823 0.215740x10" • 6 0.415225 188.607 36.1283 29.4749 ana 0.0417002 0.216120x10" -4 0.414477 186.786 36.1304 29.4772 iso 0.0417370 0.215174x10" •4 0.414477 186.965 36.1312 29.4758 siml 0.0418227 0.215352x10" -4 0.415235 185.586 36.1245 29.2810 sim2 0.0418023 0.214878x10" -4 0.415235 186.973 36.1259 29.4688 60 sim3 0.0417662 0.215733x10" -4 0.415235 186.794 36.1275 29.4741 sim4 0.0417665 0.215737x10" -4 0.415235 186.790 36.1280 29.4746 sim5 0.0417665 0.215737x10" •4 0.415235 186.790 36.1281 29.4747 sim6 0.0417665 0.215737x10" -4 0.415235 186.790 36.1281 29.4748 ana 0.100575 0.00215824 0.415564 160.987 36.1098 29.4672 iso 0.100431 0.00214979 0.415553 160.933 36.1141 29.4676 siml 0.100645 0.00215072 0.416317 160.109 36.1049 29.2716 sim2 0.100512 0.00214678 0.416310 160.938 36.1085 29.4605 600 sim3 0.100645 0.00215449 0.416320 161.067 36.1070 29.4642 sim4 0.100658 0.00215453 0.416320 161.005 36.1074 29.4647 sim5 0.100643 0.00215453 0.416320 161.003 36.1075 29.4648 sim6 0.100666 0.00215453 0.416320 161.005 36.1076 29.4648 ana 0.683376 0.190695 0.512251 141.923 34.2940 28.5883 iso 0.682959 0.191040 0.512392 141.926 34.3098 28.5986 siml 0.683443 0.191763 0.513133 141.514 34.3483 28.4171 sim2 0.682408 0.190726 0.513005 141.928 34.3066 28.5928 6000 sim3 0.684061 0.190419 0.512891 141.999 34.2956 28.5911 sim4 0.682981 0.190437 0.512876 141.945 34.2972 28.5885 sim5 0.682998 0.190439 0.512877 141.939 34.2974 28.5887 sim6 0.683092 0.190439 0.512877 141.939 34.2975 28.5887 ana 4.55387 1.70847 1.64196 110.103 21.5995 21.6613 iso 4.54685 1.70634 1.64644 110.053 21.5830 21.6629 siml 4.56176 1.71946 1.63998 109.402 21.4431 21.5415 sim2 4.55291 1.70905 1.64730 110.056 21.5819 21.6573 60000 sim3 4.57740 1.71047 1.64424 110.144 21.6141 21.6731 sim4 4.55865 1.71048 1.64354 110.105 21.6007 21.6620 sim5 4.55860 1.71067 1.64278 110.109 21.5999 21.6622 sim6 4.55913 1.71055 1.64367 110.109 21.6004 21.6619 ana 13.9902 5.11638 5.11640 102.208 18.7503 18.7503 iso 13.9102 5.08706 5.08729 102.188 18.7471 18.7467 siml 13.8688 5.08011 5.07796 101.775 18.7058 18.7052 sim2 13.9682 5.10549 5.10596 102.169 18.7377 18.7375 600000 sim3 14.0963 5.12703 5.13941 102.214 18.7554 18.7562 sim4 14.0001 5.12097 5.12041 102.203 18.7487 18.7487 sim5 14.0028 5.12067 5.12070 102.204 18.7476 18.7476 sim6 14.0079 5.12252 5.12254 102.204 18.7475 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 wedge-shaped 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°. Chapter 3. Element Types and Shape Functions 52 0.020 0.015 d o.oio-i I 0.005 o.ooo- •= -0.005 • = ana V = iso siml A — sim2 + = sim3 X sim4 8 (mm) 0.020 0.015 o.oio-0.005-12 -0.005 • — ana V — iso — siml A sim2 + = sim3 X sim4 0.000\ [nr-JW" ~* MH •'«ITM"igrJ r (mm) 0.020 0.015 I 0.020 -0.005 0.015 0.010 0.005 -0.005 • — ana 7 — iso -0 — siml A = sim2 + = aim 3 X = sim4 -c 10.0 a 5.0-0.0 -5.0 -10.0 -15.0 jag-"* -MP' <^ a = ana ^ = iso -» = siml A = sim2 + = sim3 x = sim4 18 19 i 20 r (mm) 21 22 Figure 3.11: Current density distribution in the SC coaxial cable at 6kHz Chapter 3. Element Types and Shape Functions 53 Figure 3.12: Current density distribution in the SC coaxial cable at 60kHz Chapter 3. Element Types and Shape Functions 54 46 80 0 (degree) 80 £2. d H 20 sim3 -o = 6 Hz = 60 Hz -•» = 600 Hz = 6 kHz -+ = 60 kHz -« = 600 kHz 30 46 ao 0 (degree) 20 a E o 9 16 46 0 (degree) sim3 —a = 6 Hz —v = 60 Hz -~» = 600 Hz - -A = 6 kHz •-+ = 60 kHz —« = 600 kHz 0 (degree) Figure 3.13: Maximum errors in [R] and [L] at different 0 Chapter 3. Element Types and Shape Functions 55 3.6.2 Computation efficiency: CPU 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 the final equations depicted in Fig. 3.7, type of algorithm for solving the final equations, 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 isoparametric 2nd order simplex (a) meshes with 2% error (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 pur-pose 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 parame-ters 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 Chapter 3. Element Types and Shape Functions 56 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 with floating point accelerator. As a variable bandwidth Table 3.7: Storage and other parameters for different elements element type e (degree) M N BW matrix dimension storage (bytes) Max errors (%) R L iso 90 32 89 17 1088x1 83900 0.85 0.06 siml 22.5 464 225 20 3666x1 227084 1.87 1.59 sim2 18 300 581 119 28474x1 1372376 1.65 0.52 sim3 18 220 961 206 70636x1 3838008 1.81 0.25 sim4 22.5 176 1377 319 151385x1 7387372 — — Table 3.8: CPU time requirements for different elements element type matrix formation matrix factorization solution loss-energy method others total CPU time ISO 2.4 s (26.0%) 1.8 s (19.5%) 0.5 s ( 5.4%) 3.2 s (34.3%) 1.4 s (14.8%) 9.3 s siml 2.3 s (12.5%) 6.7 s (36.8%) 1.6 s ( 8.5%) 3.5 s (18.9%) 4.3 s (23.4%) 18.4 s sim2 4.6 s ( 2.9%) 130.4 s (82.1%) 10.4 s ( 6.5%) 7.1 s ( 4.5%) 6.4 s ( 4.0%) 158.9 s sim3 11.0 s ( 1.9%) 506.7 s (88.9%) 26.7 s ( 4.7%) 11.9 s ( 2.1%) 13.9 s ( 2.4%) 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 Chapter 3. Element Types and Shape Functions 57 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 matrix dimension storage (bytes) CPU time factorization solution iso 49x89 136268 2.8 s 0.7 s siml 58x225 377228 10.7 s 2.2 s sim2 355x581 4216872 403.9 s 28.9 s sim3 616x961 12179448 1703.7 s 103.1 s sim4 955x1377 26005772 — — 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 e (degree) M N BW matrix dimension storage (bytes) Max errors (%) R L ISO 180 16 45 9 303x1 33928 10.10 2.09 siml 60 174 85 9 611x1 55604 9.75 8.06 sim2 60 90 175 29 3158x1 170384 14.06 6.45 sim3 60 66 289 61 9628x1 454104 14.30 6.14 sim4 60 66 517 105 27845x1 1230332 14.14 5.81 sim5 60 54 661 161 50056x1 2214940 13.70 5.59 sim6 60 54 955 229 99610x1 4335336 13.43 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 sim-plex 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 matrix matrix loss-energy total type formation factorization solution method others CPU time ISO 1.2 s 0.5 s 0.2 s 1.7 s 0.8 s 4.4 s (28.2%) (10.6%) ( 4.6%) (38.7%) (17.8%) siml 1.0 s 0.9 s 0.4 s 1.3 s 1.7 s 5.3 s (18.3%) (16.7%) ( 6.8%) (25.5%) (32.7%) sim2 1.1 s 6.2 s 1.3 s 2.2 s 2.0 s 12.8 s ( 8.9%) (48.2%) (10.3%) (17.3%) (15.2%) sim3 2.1 s 32.7 s 3.7 s 3.7 s 3.3 s 45.5 s ( 4.5%) (71.9%) ( 8.2%) ( 8.2%) ( 7.2%) sim4 4.8 s 154.0 s 10.9 s 7.9 s 7.6 s 185.2 s ( 2.6%) (83.2%) ( 5.9%) ( 4.3%) ( 4.1%) sim5 8.6 s 367.5 s 18.4 s 11.1 s 16.1 s 421.7 s ( 2.0%) (87.1%) ( 4.4%) ( 2.6%) ( 3.8%) sim6 17.2 s 1017.2 s 39.7 s 19.9 s 36.8 s 1130.8 s ( 1.5%) (90.0%) ( 3.5%) ( 1.8%) ( 3.3%) Table 3.12: Storage and CPU time requirements for pivoting element type matrix dimension storage (bytes) CPU time factorization solution iso 25x45 47080 0.6 s 0.3 s siml 25x85 79828 1.1 s 0.5 s sim2 85x175 357856 11.8 s 2.2 s sim3 181x289 1137000 66.9 s 6.6 s sim4 313x517 3373948 333.0 s 23.1 s sim5 481x661 6501100 897.6 s 44.1 s sim6 685x955 13208376 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 loss-energy 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 the final symmetric 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 un-derneath 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 distribu-tion 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 FEM 61 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 a field truncation boundary, and FEM with such a boundary shall be called "conventional FEM" here. With the help of the penetration depth, the field truncation 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 the field truncation 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 to find suitable 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 the field solution of a current filament. 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 the field solution of the equivalent current filament. The proposed technique creates a much smaller fixed solution 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 the field in the earth can be solved theoretically. Therefore, such a cable is used first to study the division schemes and the locations of the field truncation boundary in the earth for the conventional FEM. The impedances of the cable are calculated by the FEM with the field truncation 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 the field truncation 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,/5e and re/Se. The parameters r-fc and re are the boundary radius and inner earth radius, respectively. For re=24mm, accurate results can be obtained if rt,/8e < 0.2. 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 r0 > 12Se is required for the conventional FEM. Comparisons 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 e if the earth penetration is large. 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 the field solution of an equivalent current filament by applying the perturbation concept. For deeply buried SC coaxial cables, the earth can be assumed to occupy the whole space, and the corresponding fields become axisymmetrical. This leads Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 63 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 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 64 x-y axes as shown in the figure, the E field distributions are [2] [10] jwu.J f°» e-va-fcv'«*+i/p2 ^ = _ j a ^ i f°°_e * Cos(sa)<fa y>0 (4.1) 7T 7o a + ^ ^ a 2 + 1 / p 2 Ee = 2TT - K 0 — - K 0 — + / . cos{xa)da \ \PeJ \pj Jo e±a + ^ a 2 + i / p 2 y < 0 (4.2) where = v/x2 + (y + /i)2 (4.3) £' = y/x* + (y-hy (4.4) Pe = (4.5) Ea and J?e are the E fields in the air and in the earth, respectively, and their detailed derivations are given in Appendix B. h is the burial depth of the cable. I is the current in the filament. pe is the complex penetration depth in the earth. Ko is the zero order second kind modified Bessel function. fia is the permeability of the air. /xe and pe are, 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 the field exists. 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 fie = / i 0 and typical value of pe = 10012m for the earth, the real penetration depth in the earth Se is 5033 m at 1 Hz and 5.03 m at 1 MHz. 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 FEM 65 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, Ee(xp, yp), the E value at point P on the inner surface of the earth as shown in Fig. 4.1(a), can be calculated approximately from (4.2). The surface impedance at point P is defined as Ee(xp,yP) ZE = I = w« ^ ^ + (yp + h)2 j K u ^ yte + (yp ~ h y y + / . = cos(xp ct)da 1 (4.6) Jo + + / Strictly speaking, ZE would be different at different locations on the inner surface of the earth. The differences, however, are very small, and can be ignored. If xp = re and yp = —h, (4.6) becomes Z.= , ( K o ( - ) - K 0 ( £ ± ^ \ + f c o s ( r « a ) J 2^ r \ \peJ \ pe ) Jo ^ a + yja* + 1/pl J (4.7) re is the inner surface radius of the earth. ZE is also called "earth return impedance." Equation (4.7) was first derived by Pollaczek, and shall be called "Pollaczek's formula." Once ZE is known, [Z] of a buried SC coaxial cable can be calculated with (2.68), except that ZEQK in (2.72) has to be modified into ZEQK = ZBK + ZDK + ZE (4.8) ZE in (4.7) can be found by numerical integration. It should be noted from (2.69) and (2.70) that every element in matrix (2.68) has ZEQK as one of its components. Therefore, ZE must be added to every element in [Z] of a SC coaxial cable. Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 66 If a SC coaxial cable is buried deeply in the earth (h ^ > Se), the earth can be assumed 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 TBh —• oo and = r e. The resulting formula is Z e ~ 27rre M r ^ ) [™> re is the inner surface radius of the earth. The above equation is to be used in (4.8) 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 re does not exist. Also, it is not clear how the tunnel structure will affect the impedance calculation with respect to the insulation within the tunnel. Similarly, Ze for shallowly buried SC coaxial cables 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. Earth Region Reduction Technique for [Z] Calculation with FEM 67 4.3 [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, a field truncation boundary located at r& shown in Fig. 4.2(b) has to be established, which should enclose the most significant part of the field distribution. 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 68 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 - JsK+1 /<TK+I) (i = 1,2,..., K) (4.10) where Js{ and JsK+1 are the source current densities in the ith and the reference con-ductors, respectively. 0% 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 Jjt^ds]l{IR,IRi+IIiIIi) (z,j = l,2,. . . , /0 (4.11) L*i = { E R E ( / 5 C W(V 5 ) } / ( 7 ^+W/>) (*,j = l,2,...,2iQ (4.12) where 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 e can be calculated analytically. If a loop current of 1+j'O A 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 Ie from the analytical formula is shown in Fig. 4.3(a), as a function of rj, at 6kHz. Se is 64.97m at this frequency. Similar curves can be obtained for other frequencies. It can be seen that IeR wl A and Iei «0 at 5Se. This means that the most significant part of the field is enclosed by a boundary at 5Se for the deeply buried SC coaxial cable. FEM is applied to calculate [Z] of the system in Fig. 4.2 for different r 0. Isoparametric elements are used in the calculation, and the solution region is reduced to a wedged region Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 69 1 2 3 4 6 1 2 3 4 9 rb/6e rb/6e (a) enclosed earth current by r& at 6kHz (b) max errors in [R] and [L] 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 10n, 10n+?, 10 n +5, 10n + 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=35e. For the same system, pe is varied from lOOOfhn to O.Olfim, and the results show that 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£e can be used as the 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 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 70 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 r0. 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"+1. These patterns are listed in Tab.4.1. The first earth division radius is re and the last one is r0=3Se. Between r e and r&, the earth is divided by division radii listed in Tab.4.1. Table 4.1: Earth division radius patterns for isoparametric elements in the decade from 10" to 10"+1 pattern division radii 1 10" 10"+i 10"+* 10"+1 2 10" 10"+2 10"+1 3 10" i x 10"+1 \ x 10"+1 10n + 1 4 10" \ x 10"+1 10n + 1 5 10" 10"+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) division radii (m) 6 0.024 0.0464 0.1 0.215 0.464 1 2.15 4.64 10 21.5 46.4 100 215 464 1000 2150 4640 6164 60 0.024 0.0464 0.1 0.215 0.464 1 2.15 4.64 10 21.5 46.4 100 215 464 1000 1949 600 0.024 0.0464 0.1 0.215 0.464 1 2.15 4.64 10 21.5 46.4 100 215 464 616 6000 0.024 0.0464 0.1 0.215 0.464 1 2.15 4.64 10 21.5 46.4 100 195 60000 0.024 0.0464 0.1 0.215 0.464 1 2.15 4.64 10 21.5 46.4 61.6 600000 0.024 0.0464 0.1 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 /(Hz) element R (ft/km) L (mH/km) Ru R\2 R22 L\\ L12 L22 6 ana iso 0.0447332 0.0059220 0.420388 0.0447221 0.0059109 0.420377 2.41400 2.26152 2.25486 2.41099 2.25853 2.25188 60 ana iso 0.100918 0.0592392 0.473695 0.100821 0.0591056 0.473561 2.18191 2.03126 2.02461 2.17968 2.02884 2.02219 600 ana iso 0.692750 0.594334 1.00774 0.691494 0.593213 1.00662 1.92586 1.80098 1.79434 1.92364 1.79882 1.79218 6000 ana iso 6.60507 6.11239 6.43395 6.59130 6.09938 6.42073 1.67654 1.56891 1.56320 1.67495 1.56733 1.56162 60000 ana iso 63.7665 60.9211 60.8546 63.6482 60.8077 60.7478 1.41446 1.32596 1.32602 1.41308 1.32461 1.32469 600000 ana iso 605.816 596.942 596.942 604.401 595.578 595.578 1.17633 1.09287 1.09287 1.17555 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. re is varied among the values: 50mm, 100mm, 250mm, 500mm, and 1000mm with />e=100flm, while the other parameters remain the same. When r e is larger than 50mm, additional 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 FEM 72 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 Z5e away from the cables 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 10n, 10 n + 3 ,10 n + § , 10n + 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 the field distribution 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 the field inside 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 con-ventional 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 gen-eral 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 AH = - ~ in air (4.13) Eb 1 Aj = — : h Js in earth (4.14) ]U3(Te where Js and <re are the source current density and conductivity in the earth, respectively. 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 [AB] = JM + M ^ ± i (4.15) where [ly] is a vector with NB elements. If the boundary nodes are numbered such that the first NBB nodes are the boundary nodes in the earth, and the remaining NB—NBB nodes are in the air, the first NBB elements in [ly] will be unity and the rest of the elements will be zero. [AB] becomes a partially known vector. [IEA and [IE2] m (3.38) become fe] = -ID.] ( -J f f 1 + = [W - [ W « r „ (4-16) [/*] = MGc][FB]T ( - ^ + M ^ ± l ) = [IEI] + [SCv]JsK„ (4.17) Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 76 where [J*] = [DB][EB]/U») VBA = -[GC}[FB]T[EB) [4.10) [FV] = [Db][1V]/UV<TK+I) [Scv] = [Gc][FB]T[lv]/<rK+i [Fv] and [Scv] are vectors. Putting [J^ ,] and [IEJ] in (4.16) and (4.17) into (3.37) and moving JsK+1 terms to the left side produces [D] -[FA] ' " [Au] ' [FH] [ScA] . . [J] + [IEA] . (4.19) [F] becomes [FA] if [Fv] is subtracted from its last column, and [Sc] becomes [ScA] if [ScY] is subtracted from its last column. [I] always corresponds to a loop current. For a loop current Ij between conductor j and the earth, the jth and (iif+l)th elements in [I] will be Ij and J e p , respectively. J e p is the partial earth return current enclosed by the 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 JsK+1 in the earth, the corresponding boundary nodes are still treated as having fixed boundary values. They do not appear in the final equations of the solution as shown by (4.19). Their values will be evaluated by (4.15) after solving for JsK+1-For a deeply buried SC cable, the field solutions of both the original cable and the buried current filament have 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 current filament is drawn in Fig. 4.5(b). Ec(r) and Ep(r) in Fig. 4.5 are the E fields associated with the cable and the filament, respectively. They Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 77 tf- - earth it A / ^ S C cable / 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 Ec(r) = jW(leI K 0 {r/pe) Mr) = - ^ K 0 ( r / P e ) 2TT (r e/ P e)Ki(r e/p e) r > 0 r >re 2TT (4.20) (4.21) I is the loop current between the cable or the filament and the earth. The partial earth return currents between radii r e and r for both cases in Fig. 4.5 are (TIVe) \ {T'IP*)t r > r e ) W . ) = M ^ ) - ^ ) = - ( r e / p e ) K J i ( r e / p e ) ( i - ^ ^ (4.22) (4.23) where IF(r) = Iepc(r)\r^Q = -1(1 - (r/p^Mr/p,)) (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 ( ? ) <4-25> Pe \PeJ Ec(r) and Jepc(r) can be respectively related to Ep(r) and Ieppi?) as £7c(r) = E^r)/^ r >re (4.26) Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 78 Iepc(r) = CpIepF(r) r > r e (4.27) The values of Cp at different |r«/pe| are plotted in Fig. 4.6. From Fig. 4.6 it can be seen that Cp is close to 1.0 when the earth penetration is large compared with r e, since l-r 0.989-^ 0.998-0.997-0.996-0.00 0.02 0.04 0.06 0.08 0.10 0.00 0.02 0.04 0.06 0.08 0.10 K/P.\ \r./p.\ Figure 4.6: Perturbation coefficient Cp as a function of |r e/p 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 |re/pe|=0.02, /=lMHz, />e=100£2m, and /ze=/io> |pe|=3.56m and re=71.2mm. By using the exact formulas in (4.20) and (4.22) to find [EB] in (4.18) and Iep for [I] in (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 Iep for rfr=0.24m, lm, 2.5m, 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& Se. When the frequency is high enough such that Se becomes comparable to r e or even smaller than r 0, the final matrix 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 rt are tested for the above case with the r& listed in Tab.4.4, and the numerical results show that the errors are less than 1% in [R] for < S€ and in [L] for r& < 2Se. pe is varied among lOOOftm, lOflm, lflm, O.lfim, and O.Olfim with re=24mm, while r e is varied among 50mm, 100mm, 250mm, 500mm, 1000mm, and 2000mm with />e=100fim. Table 4.4: [Z] of the deeply buried cable found from a reduced earth region / n R (n/km) L (mH/km) (Hz) (m) Rn Rl2 #22 Irn Xl2 ana 0.0447332 0.00592198 0.420388 2.41400 2.26152 2.25486 0.24 0.0447332 0.00592198 0.420388 2.41323 2.26079 2.25413 1.00 0.0447332 0.00592200 0.420388 2.41274 2.26028 2.25362 6 2.50 0.0447332 0.00592198 0.420388 2.41247 2.26001 2.25336 10.00 0.0447332 0.00592198 0.420388 2.41191 2.25945 2.25280 25.00 0.0447332 0.00592198 0.420388 2.41164 2.25918 2.25253 ana 0.100918 0.0592392 0.473695 2.18191 2.03126 2.02461 0.24 0.100954 0.0592391 0.473695 2.18138 2.03055 2.02390 1.00 0.100955 0.0592392 0.473695 2.18086 2.03002 2.02337 60 2.50 0.100955 0.0592393 0.473695 2.18059 2.02975 2.02310 10.00 0.100955 0.0592391 0.473695 2.18003 2.02920 2.02254 25.00 0.100955 0.0592392 0.473695 2.17976 2.02892 2.02227 ana 0.692750 0.594334 1.00774 1.92586 1.80098 1.79434 0.24 0.692605 0.594326 1.00773 1.92508 1.80027 1.79362 1.00 0.692607 0.594326 1.00773 1.92457 1.79975 1.79310 600 2.50 0.692607 0.594326 1.00773 1.92430 1.79948 1.79283 10.00 0.692606 0.594325 1.00773 1.92374 1.79892 1.79228 25.00 0.692577 0.594295 1.00770 1.92347 1.79865 1.79200 ana 6.60507 6.11239 6.43395 1.67654 1.56891 1.56320 0.24 6.60474 6.11280 6.43415 1.67583 1.56821 1.56250 1.00 6.60466 6.11274 6.43409 1.67530 1.56768 1.56197 6000 2.50 6.60467 6.11275 6.43410 1.67503 1.56741 1.56170 10.00 6.60393 6.11204 6.43340 1.67448 1.56686 1.56115 25.00 6.60439 6.11248 6.43383 1.67421 1.56659 1.56088 ana 63.7665 60.9211 60.8546 1.41446 1.32596 1.32602 0.24 63.7597 60.9191 60.8592 1.41369 1.32523 1.32531 1.00 63.7593 60.9189 60.8590 1.41317 1.32470 1.32478 60000 2.50 63.7588 60.9184 60.8585 1.41290 1.32443 1.32451 10.00 63.7098 60.8736 60.8137 1.41235 1.32388 1.32396 25.00 63.8667 61.0229 60.9630 1.41131 1.32287 1.32295 ana 605.816 596.942 596.942 1.17633 1.09287 1.09287 0.24 605.736 596.914 596.914 1.17559 1.09215 1.09215 1.00 605.711 596.892 596.892 1.17507 1.09163 1.09163 600000 2.50 605.635 596.826 596.826 1.17480 1.09136 1.09136 10.00 598.441 590.205 590.205 1.17388 1.09047 1.09047 25.00 586.913 579.558 579.558 1.18241 1.09834 1.09834 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 80 The earth is meshed with the pattern 10n, 10n+*, 10n+*, 10n + 1. The numerical results 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 approx-imate Eb and Iep from the filament field solutions for different r 6 are given in Tab.4.5. Iep is calculated by numerical integration based on Ep(r). The FEM solution meshes are used for the numerical integration as well. With re=24mm, the results show a similar pattern as with the exact Eb and Iep: accurate impedances are obtained for the r e listed in Tab.4.5 if 8e is much larger than r 6, and erroneous impedances are obtained if Se is comparable with rb. For the rb listed in Tab.4.5, if the earth resistivity is varied among lOOOfim, lOfim, lfim, O.lfim, and O.Olfhn with the same r e, the errors are less than 1% in [R] for r 0 < 0.2Se and in [L] for rb < 0.5£e. For the same Se, the larger the r e, the larger the ratio re/Se, and the larger the errors in Eb and Iep given by the filament formulas. Therefore, for a large r e, the ratio re/Se has to be less than a certain value in order to achieve accurate results with the proposed technique. By varying re among 50mm, 100mm, 250mm, 500mm, 1000mm, and 2000mm with />=100ftm, the results show that the errors are less than 1% in [R] for re/Se < 0.012, and in [L] for re/Se < 0.055. 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 8e for small r e, and by the relationship between r e and 8e for large r e . For the cable with />=100fhn, the ratio re/8e becomes more influential in determining the frequency range if r e > 250mm. Table 4.5: [Z] found with the proposed technique based on Ep f n R (ft/km) L (mH/km) (Hz) (m) Rii R\2 R22 Lu L12 L22 ana 0.0447332 0.00592198 0.420388 2.41400 2.26152 2.25486 0.24 0.0447332 0.00592197 0.420388 2.41322 2.26079 2.25413 1.00 0.0447332 0.00592198 0.420388 2.41274 2.26028 2.25362 6 2.50 0.0447332 0.00592198 0.420388 2.41247 2.26001 2.25336 10.00 0.0447332 0.00592198 0.420388 2.41191 2.25945 2.25280 25.00 0.0447332 0.00592198 0.420388 2.41164 2.25918 2.25253 ana 0.100918 0.0592392 0.473695 2.18191 2.03126 2.02461 0.24 0.100954 0.0592391 0.473695 2.18138 2.03055 2.02390 1.00 0.100955 0.0592391 0.473695 2.18086 2.03002 2.02337 60 2.50 0.100955 0.0592391 0.473695 2.18059 2.02975 2.02310 10.00 0.100955 0.0592391 0.473695 2.18003 2.02920 2.02254 25.00 0.100955 0.0592391 0.473695 2.17976 2.02893 2.02227 ana 0.692750 0.594334 1.00774 1.92586 1.80098 1.79434 0.24 0.692605 0.594326 1.00773 1.92508 1.80027 1.79362 1.00 0.692607 0.594326 1.00773 1.92457 1.79975 1.79310 600 2.50 0.692607 0.594326 1.00773 1.92430 1.79948 1.79283 10.00 0.692606 0.594324 1.00773 1.92374 1.79892 1.79228 25.00 0.692575 0.594292 1.00770 1.92347 1.79866 1.79201 ana 6.60507 6.11239 6.43395 1.67654 1.56891 1.56320 0.24 6.60480 6.11285 6.43420 1.67583 1.56821 1.56250 1.00 6.60470 6.11278 6.43414 1.67530 1.56768 1.56197 6000 2.50 6.60472 6.11280 6.43415 1.67503 1.56741 1.56170 10.00 6.60399 6.11210 6.43345 1.67448 1.56686 1.56115 25.00 6.60443 6.11252 6.43387 1.67421 1.56659 1.56088 ana 63.7665 60.9211 60.8546 1.41446 1.32596 1.32602 0.24 63.7646 60.9237 60.8638 1.41369 1.32522 1.32530 1.00 63.7643 60.9235 60.8637 1.41317 1.32470 1.32478 60000 2.50 63.7633 60.9226 60.8627 1.41290 1.32443 1.32451 10.00 63.7145 60.8780 60.8181 1.41235 1.32388 1.32396 25.00 63.8656 61.0218 60.9617 1.41132 1.32287 1.32295 ana 605.816 596.942 596.942 1.17633 1.09287 1.09287 0.24 606.086 597.237 597.237 1.17556 1.09213 1.09213 1.00 606.059 597.214 597.214 1.17504 1.09160 1.09160 600000 2.50 605.980 597.146 597.146 1.17477 1.09133 1.09133 10.00 598.766 590.504 590.504 1.17386 1.09045 1.09045 25.00 586.921 579.565 579.566 1.18240 1.09833 1.09833 If re=500mm with p—lOOilm, the minimum 8e will be 41.7m, and the corresponding frequency is /=pe/(7r/zei52)«15kHz. Therefore, the frequency range for the proposed Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 82 technique is from 1Hz to 15kHz, and rj, is assigned a fixed value, say 2.5m, for this frequency range. For the conventional FEM, r0>36e is required for the deeply buried 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 pre-ceding 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 e. The impedances of a SC cable with an 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 r0 > 12Se will give reasonably accurate results. Based on the comparison with the conventional FEM it is shown that accurate results can be obtained with the proposed technique for different pe and r e. Good agreements are observed among the field distributions around 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 re/Se is small. For re=24mm, the maximum differences between the results from Pollaczek's formula and those from the FEM are less than 1% in R if re/Se < 0.03 and in L if rJ8e <0.095. For r„=1000mm, the maximum differences are less than 1% in R if re/Se <0.018 and in L if re/^e <0.154. The numerical results for the cable with an arbitrary structure show that Pollaczek's formula can still be used with an approximate r e. Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 83 4.5.1 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.5m and re=24mm in Fig. 4.1(a). It is also assumed that />e=100f2m and / z e = / i 0 . Once 77, is chosen, the earth is divided with the isoparametric elements in the pattern 10n, 10n+3, 10"+t, 10n + 1 discussed in Section 4.3. The mesh for rj,=2£e at 6kHz is shown in Fig. 4.7. (a) entire mesh (b) detailed mesh around the cable Figure 4.7: FEM mesh at 6kHz for rb=2Se The impedances calculated with r& being equal to 55e, 106,., and 15Se, and with Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 84 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=5Se. For [R] good agreement exists at rb=15Se. Table 4.6: [Z] of the shallowly buried SC coaxial cable from the conventional FEM / (Hz) n R (fi/km) L (mH/km) Rll Rl2 R-22 L\\ L\2 l>22 6 56, 106, 156, Pollaczek 0.0445242 0.00570694 0.420238 0.0446705 0.00585330 0.420384 0.0447053 0.00588803 0.420419 0.0447405 0.00592928 0.420395 2.50775 2.35534 2.34869 2.50897 2.35656 2.34990 2.50916 2.35675 2.35010 2.51380 2.36132 2.35467 60 56, 106, 156, Pollaczek 0.0989151 0.0571899 0.471710 0.1004902 0.0587650 0.473285 0.1008098 0.0590845 0.473605 0.1011468 0.0594682 0.473924 2.27548 2.12469 2.11803 2.27677 2.12598 2.11932 2.27695 2.12615 2.11950 2.28130 2.13064 2.12399 600 56, 106, 156, Pollaczek 0.677598 0.579357 0.992824 0.692109 0.593867 1.007335 0.695547 0.597305 1.010773 0.699820 0.601403 1.014809 2.01881 1.89399 1.88734 2.01996 1.89514 1.88849 2.02014 1.89532 1.88867 2.02392 1.89904 1.89240 6000 56, 106, 156, Pollaczek 6.59260 6.10068 6.42211 6.74617 6.25425 6.57568 6.77707 6.28515 6.60658 6.81484 6.32215 6.64371 1.76564 1.65802 1.65231 1.76670 1.65908 1.65337 1.76683 1.65921 1.65350 1.77047 1.66284 1.65713 60000 56, 106, 156, Pollaczek 67.3478 64.5070 64.4471 68.6650 65.8243 65.7643 68.9727 66.1319 66.0720 69.3549 66.5095 66.4430 1.49249 1.40402 1.40410 1.49300 1.40453 1.40461 1.49306 1.40458 1.40466 1.49589 1.40739 1.40745 600000 56, 106, 156, Pollaczek 698.632 689.808 689.808 709.490 700.667 700.667 711.493 702.670 702.670 714.613 705.740 705.740 1.22252 1.13908 1.13908 1.22201 1.13856 1.13856 1.22186 1.13842 1.13842 1.22453 1.14107 1.14107 The earth return current enclosed by the boundary is calculated from Ee in (4.2) with 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,=nSe and the [Z] calculated with r(,=(n-l)5e is denned as the nth maximum difference. It is taken from all the elements in [Z] in the frequency range from 1Hz to Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 85 rb/se rb/6e (a) enclosed earth current by r0 at 6kHz (b) max differences in [R] and [L] Figure 4.8: Earth return current and maximum differences in [Z] for different r0 1MHz. Fig. 4.8 shows that the maximum difference in [R] decreases slowly. It reaches 0.07% at r0=15Se. The maximum difference in [X] reaches 0.027% even at r0=6Se. If the nth maximum difference is defined as the one between r6=(ra-3).5e and r0=nSe and if 0.5% is chosen as the threshold value of the maximum difference in [R] for stopping the iteration, r0=12Se will be the iteration result because the maximum difference in [R] is 0.55% at 12Se and .27% at 15Se. The numerical results show that the maximum difference in [R] between 12^ and 158e maintains the value 0.27% at different pe and r e. For pe=100flm and re=24mm, the differences between the results found with the FEM having r0=12Se and those found with Pollaczek's formula are less than 1.1% in [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, >12Se will give reasonably accurate impedances for the shallowly buried SC cables. Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 86 4.5.2 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, re=24mm, />e=100fhn, and pe=po. Four different r0 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 r0 are listed in Tab.4.7. "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 Iep 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 E field distributions 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 r0=5m, and "analytical" represents the 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. Iep 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 87 Table 4.7: [Z] of the shallowly buried SC coaxial cable from the proposed technique / n R (ft/km) L (mH/km) (Hz) (m) R\\ R\2 #22 In 2/12 L22 conventional 0.0447186 0.00590137 0.420432 2.50922 2.35681 2.35015 2.5 0.0447465 0.00592928 0.420460 2.51230 2.35989 2.35323 5 0.0447465 0.00592928 0.420460 2.51207 2.35966 2.35300 6 10 0.0447465 0.00592928 0.420460 2.51172 2.35931 2.35266 25 0.0447465 0.00592928 0.420460 2.51144 2.35903 2.35238 Pollaczek 0.0447405 0.00592928 0.420395 2.51380 2.36132 2.35467 conventional 0.100950 0.0592246 0.473745 2.27700 2.12620 2.11955 2.5 0.101193 0.0594680 0.473988 2.28000 2.12921 2.12255 5 0.101193 0.0594680 0.473988 2.27977 2.12898 2.12232 60 10 0.101193 0.0594680 0.473988 2.27943 2.12863 2.12198 25 0.101193 0.0594681 0.473988 2.27915 2.12835 2.12170 Pollaczek 0.101147 0.0594682 0.473924 2.28130 2.13064 2.12399 conventional 0.696866 0.598624 1.01209 2.02019 1.89537 1.88872 2.5 0.699635 0.601394 1.01486 2.02243 1.89761 1.89096 5 0.699635 0.601394 1.01486 2.02220 1.89738 1.89073 600 10 0.699636 0.601394 1.01486 2.02185 1.89703 1.89039 25 0.699640 0.601399 1.01487 2.02157 1.89675 1.89011 Pollaczek 0.699820 0.601403 1.01481 2.02392 1.89904 1.89240 conventional 6.79056 6.29865 6.62008 1.76687 1.65925 1.65354 2.5 6.81436 6.32245 6.64388 1.76903 1.66142 1.65570 5 6.81438 6.32246 6.64389 1.76880 1.66119 1.65547 6000 10 6.81443 6.32251 6.64394 1.76846 1.66084 1.65513 25 6.81482 6.32290 6.64433 1.76818 1.66056 1.65485 Pollaczek 6.81484 6.32215 6.64371 1.77047 1.66284 1.65713 conventional 69.0895 66.2488 66.1888 1.49306 1.40459 1.40467 2.5 69.3455 66.5047 66.4448 1.49441 1.40594 1.40602 5 69.3470 66.5063 66.4463 1.49418 1.40571 1.40579 60000 10 69.3523 66.5115 66.4516 1.49383 1.40536 1.40544 25 69.3739 66.5332 66.4732 1.49350 1.40503 1.40511 Pollaczek 69.3549 66.5095 66.4430 1.49589 1.40739 1.40745 conventional 712.336 703.512 703.512 1.22179 1.13835 1.13835 2.5 714.293 705.469 705.469 1.22308 1.13964 1.13964 5 714.421 705.597 705.597 1.22284 1.13940 1.13940 600000 10 714.270 705.447 705.447 1.22235 1.13891 1.13891 25 713.826 705.003 705.003 1.22215 1.13871 1.13871 Pollaczek 714.613 705.740 705.740 1.22453 1.14107 1.14107 Table 4.8: Iep found with numerical integrations for the proposed technique / (Hz) Re(Jep) and Im(/ep) (p.u.) r(,=2.5m rj,=5m rj,=10m r t t=25m 60 Re(Je„) Im(I„) 0.583064xl0-5 0.486857xl0-4 0.233561 xlO"4 0.174161 xlO"3 0.936563xl0"4 0.603340xl0-3 0.589085xl0-3 0.302400xl0-2 6000 Re(/ep) Im(I«„) 0.598188 xlO"3 0.314763xl0"2 0.241264xl0"2 0.105046X10-1 0.973100xl0"2 0.324944 x 10"1 0.598538X10"1 0.125884 600000 Im(J„) 0.641595xl0-x 0.130791 0.233091 0.310736 0.646749 0.434091 1.005062 0.895586X10-1 JR (A/km2) (A/km2) J/.=250.20+17.26-(t-l) conventional Horizontal distance (m) analytical Figure 4.10: J distributions in the earth at 6kHz from the three approaches JA=250.95+17.25 (t-l) proposed (r0=5m) Chapter 4. Earth Region Reduction Technique for [Z] Calculation vrith FEM 89 Table 4.9: Storage and other parameters for the proposed technique (rj,=5m) / (Hz) M N BW matrix dimension storage (bytes) 60 283 88 36 4351x1 302064 6k 343 106 36 5356x1 365976 600k 403 124 36 6361x1 429888 Table 4.10: CPU time requirements for the proposed technique (r4=5m) / (Hz) matrix formation matrix factorization solution others total total-/,,, 60 6.8 s ( 8.6%) 9.9 s (12.4%) 8.0 s (10.1%) 4.2 s ( 5.3%) 50.4 s (63.5%) 79.3 s 28.9 s 6k 8.7 s ( 9.9%) 11.8 s (13.4%) 7.7 s ( 8.7%) 5.1 s ( 5.7%) 55.1 s (62.3%) 88.4 s 33.3 s 600k 10.1 s (11.3%) 13.8 s (15.4%) 8.1 s ( 9.1%) 5.8 s ( 6.5%) 51.8 s (57.8%) 89.6 s 37.8 s For the conventional FEM with r0=12Se) the storage and CPU time requirements for 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,=126e) f (Hz) M N BW matrix dimension storage (bytes) 60 463 142 32 9196 x 1 507600 6k 463 142 32 9196x 1 509328 600k 463 142 33 9259x 1 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 pe, pt% and r e, as well as the locations of the equivalent current filaments do not change, Iep only needs to be calculated once. By deducting the time for calculating Iep from the total CPU time, the proposed technique will require less CPU time than the conventional FEM. Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 90 Table 4.12: CPU time requirements for the conventional FEM (r0=12Se) f (Hz) matrix formation matrix factorization solution others total 60 11.3 s (25.8%) 21.9 s (50.0%) 4.1 s ( 9.4%) 6.5 s (14.7%) 43.8 s 6k 11.6 s (27.0%) 20.7 s (48.4%) 4.0 s ( 9.3%) 6.5 s (15.2%) 42.8 s 600k 11.9 s (27.2%) 21.3 s (48.8%) 3.9 s ( 9.0%) 6.5 s (15.0%) 43.6 s If pe is varied with r e remaining at 24mm, the maximum differences between the re-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 pe Pe (fim) rj=2.5m r»=5m r6=10m rj=25m in R in L in R in L in R in L in R in L 1000 1.05% 0.15% 1.05% 0.14% 1.05% 0.13% 1.05% 0.12% 100 1.05% 0.15% 1.05% 0.14% 1.05% 0.12% 1.05% 0.11% 10 1.05% 0.14% 1.05% 0.13% 1.05% 0.12% 1.05% 0.10% 1 1.02% 0.14% 1.02% 0.13% 1.02% 0.11% 1.02% 0.09% 0.1 1.02% 0.13% 1.02% 0.12% 1.02% 0.10% 1.02% 0.08% If re is varied with pe remaining at lOOflm, the corresponding maximum differences are listed in Tab.4.14 Table 4.14: Maximum differences in [Z] with the proposed technique at different r e r, (mm) rj=2.5m rj,=5m r6=10m r»=25m in R in L in R in L in R in L in R in L 50 1.05% 0.15% 1.05% 0.14% 1.05% 0.12% 1.05% 0.11% 100 1.05% 0.15% 1.05% 0.14% 1.05% 0.12% 1.05% 0.11% 250 1.52% 0.15% 1.23% 0.14% 1.13% 0.13% 1.05% 0.11% 500 4.11% 0.62% 2.89% 0.53% 1.91% 0.45% 1.23% 0.15% 1000 — — 7.99% 1.72% 4.56% 1.42% 2.20% 0.48% From the results presented in this subsection it can be seen that for small re/8e 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 Iep is calculated only once, it will require less time. 4.5.3 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 dif-ferent 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 pe and r e. The maximum differences between 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 pe Pe (Clm) r,=24mm r«=100mm re=250m re=1000mm in R in L in R in L in R in L in R in L 1000 0.48% 0.22% 0.48% 0.22% 0.64% 0.24% 4.35% 0.26% 100 0.48% 0.24% 0.97% 0.22% 2.80% 0.23% 13.35% 2.20% 10 0.48% 0.25% 3.28% 0.25% 8.85% 1.57% 16.49% 8.50% 1 2.56% 0.24% 12.36% 2.91% 20.29% 8.22% 172.49% 9.71% 0.1 9.66% 3.29% 20.34% 15.98% 47.96% 15.47% 14394.04% 9.66% 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 pe and re are less than 20%. If pe, fie, and re are fixed, there is generally a frequency beyond which the differences become larger than a specified tolerance value. This frequency shall be called "threshold frequency." For a 1% tolerance, the threshold frequencies for different pe and re are listed in Tab.4.16. By converting the threshold frequencies together with the corresponding pe, fie, and r e into a parameter re/8e, it is found out that for a given r e, this parameter re/6e Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 92 is almost constant for different pe. With re=24mm, re/Se is 0.0213 in [R] and 0.0674 in [L]. This means that if re/Se <0.0213, the differences between Pollaczek's formula and FEM are less than 1%. re/8e related to each r e is listed at the bottom of Tab.4.16. Table 4.16: Threshold / with Pollaczek's formula for maximum differences <1% Pe rs=24mm re=100mm r,=250m r,=1000mm in R in L in R in L in R in L in R in L 1000 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 60kHz > 1MHz 100 > 1MHz > 1MHz > 1MHz > 1MHz 100kHz > 1MHz 6kHz 400kHz 10 > 1MHz > 1MHz 100kHz > 1MHz 10kHz 400kHz 600Hz 40kHz 1 200kHz > 1MHz 10kHz 200kHz 1kHz 40kHz 60Hz 4kHz 0.1 20kHz 200kHz 1kHz 20kHz 100Hz 4kHz 6Hz 400Hz r./6. 0.02133 0.06744 0.01987 0.08886 0.01571 0.09935 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) r,=24rrun re=100mm r,=250m r«,=1000mm in R in L in R in L in R in L in R in L 1000 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 100 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 400kHz > 1MHz 10 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 40kHz > 1MHz 1 > 1MHz > 1MHz 600kHz > 1MHz 100kHz > 1MHz 4kHz > 1MHz 0.1 > 1MHz > 1MHz 60kHz 400kHz 10kHz 100kHz 400Hz > 1MHz rt/6t — — 0.15391 0.39738 0.15708 0.49673 0.12566 — Table 4.18: Threshold / with Pollaczek's formula for maximum differences<30% Pe (fim) re=24mm re=100mm re=250m re=1000mm in R in L in R in L in R in L in R in L 1000 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 100 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 10 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 1 > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz > 1MHz 400kHz > 1MHz 0.1 > 1MHz > 1MHz > 1MHz > 1MHz 800kHz > 1MHz 40kHz > 1MHz re/6e — — — — 1.40496 — 1.25664 — 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 Chapter 4. Earth Region Reduction Technique for [Z] Calculation with FEM 93 results for pe in the range of lOOOfim to lfim, and for r e in the range of 24mm to 250mm, 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, /i2=0.25m, io=0.5m, pe—lflm, and p,e=p,0. The impedances of the cable are calculated with the conventional FEM and with the proposed technique. As there is no definite r e here, Pollaczek's formula is used with two approximate r e of 250mm and 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 (rfc=5m) Figure 4.11: The layout of a tunnel installed cable and the FEM mesh at 600kHz The impedances found with three approaches are listed in Tab.4.19. "Pollaczekl" represents the results from Pollaczek's formula with re=250mm, and "Pollaczek2" rep-resents the results with re=500mm. Tab.4.19 indicates that for small re/8e, Pollaczek's formula gives results which are almost independent of r e. The maximum differences be-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 94 Table 4.19: [Z] of the tunnel installed cable from three approaches / (Hz) approach R (ft/km) L (mH/km) #11 #12 #22 i l l 2/12 2/22 6 conventional proposed Pollaczekl Pollaczek2 0.0447930 0.00597572 0.420506 0.0448095 0.00599228 0.420523 0.0448038 0.00599263 0.420459 0.0448038 0.00599255 0.420458 2.04825 1.89584 1.88918 2.05020 1.89779 1.89113 2.05154 1.89906 1.89241 2.05154 1.89906 1.89241 60 conventional proposed Pollaczekl Pollaczek2 0.102843 0.0611180 0.475638 0.103030 0.0613047 0.475825 0.103012 0.0613333 0.475789 0.103005 0.0613260 0.475781 1.81237 1.66158 1.65492 1.81420 1.66341 1.65675 1.81533 1.66467 1.65802 1.81533 1.66467 1.65802 600 conventional proposed Pollaczekl Pollaczek2 0.741197 0.642955 1.05642 0.745902 0.647661 1.06113 0.748330 0.649913 1.06332 0.747710 0.649293 1.06270 1.54492 1.42010 1.41345 1.54557 1.42075 1.41410 1.54678 1.42191 1.41526 1.54680 1.42193 1.41528 6000 conventional proposed Pollaczekl Pollaczek2 7.39506 6.90314 7.22457 7.50986 7.01794 7.33937 7.66821 7.17553 7.49709 7.61851 7.12583 7.44739 1.26626 1.15864 1.15293 1.26433 1.15671 1.15100 1.26433 1.15670 1.15099 1.26459 1.15696 1.15125 60000 conventional proposed Pollaczekl Pollaczek2 61.1415 58.3008 58.2408 60.5850 57.7442 57.6843 68.8430 65.9976 65.9311 65.5297 62.6843 62.6177 0.97484 0.88637 0.88645 0.96679 0.87832 0.87840 0.95450 0.86600 0.86606 0.95774 0.86923 0.86930 600000 conventional proposed Pollaczekl Pollaczek2 362.441 353.617 353.617 355.463 346.639 346.639 506.977 498.103 498.103 361.874 353.000 353.000 0.79508 0.71164 0.71164 0.79609 0.71264 0.71264 0.72688 0.64342 0.64342 0.75557 0.67211 0.67211 the conventional FEM are 30.87% in [#] and 14.25% in [L] with re=250mm, and are 14.88% in [#] and 6.95% in [L] with re=500mm. If pe is lOOfhn, the maximum differ-ences become 5.23% in [#] and 0.42% in [L] with re=250mm, and become 4.18% in [#] and 0.38% in [L] with re=500mm. These results suggest that the impedances of a cable layout of arbitrary structures can still be calculated with Pollaczek's formula by using an approximate re. 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 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 f l.=0.37006-0.01807-(»-l) / i i i=0.36736-0.01921-(«-l) 32. JA=0.83914-0.04282-(«-l) jrJi=0.84298-0.04604-(t-l) conventional 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 pene-tration 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 shal-lowly 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 0 > 3»5e. With the proposed technique, accurate results are obtained if ri,/8e is small at small r e, or if re/Se is small at large r e. For re=24mm, the errors of the proposed 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, rb > 12Se is required for the conventional FEM. Comparisons with Pollaczek's formula for the shallowly buried cables show that there are discrepancies with the con-ventional FEM when the earth penetration is small. For typical ranges of pe and r e, however, the discrepancies are reasonably small. With pe varying between lOOOJlm to lfim and with r e varying between 24mm to 250mm, the maximum discrepancies are less than 21% in [R] and less than 9% in [£]. The field solutions 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 to find impedances 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 capaci-tance calculation into a two-dimensional steady-state electric field solution 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 magnetic field can be applied directly to solve the electric field. The corresponding software can easily be adapted to handle the electric field solution. Therefore, only FEM is used in this thesis. In this chapter the general procedures for applying FEM to the steady-state electric field solution are first discussed. The energy method and the surface charge method to calculate the capacitances from the field solutions 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 = £ = i v . J = - i ^ (5-1) e <T a dt where c is permittivity, p is volume charge density, and a is conductivity. It is pointed out in [10] that, as p = poe-^/^1, any charges introduced into a conductor will dilute themselves to the surface with a time constant of e/cr. For a poor conducting earth with pe = 1000 fim, the corresponding time constant is still a very small value in the order of 10~8 s. In cable related transient studies, the smallest time step is typically in the order 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 no field inside the conductors, the solution region for the electric field will be composed of insulations bounded partly by conductor surfaces. Introducing the scalar potential <f> as E = - V 0 (5.2) the following Laplace equation can be derived eV20 = 0 (5.3) This is the governing equation describing the electric field. The boundary conditions are still Dirichlet boundary conditions on To and homogeneous Neumann boundary condi-tions 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>n<pn (5.4) n=l where 4>n is the value of <f> at FEM node n, and NT and <pn are the same as before. With the same procedure as discussed in Chapter 2, the following algebraic equation can be derived [U}[<f>] = [0] (5.5) where [d>] = [<j>i,d>2,...,d>NT]T (5.6) Umn = f eV<pm-V(pnds (m = 1,2,...,N;n = 1,2,...,NT) (5.7) JsR 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.€ifEids in SEI (5.8) n=l where <f>E* is the node value of <j> in the element, the integral Umn in 5^ becomes U*=eEi[ VvZ-Vrf'ds (5.9) JsB. Chapter 5. Admittance Calculation with Finite Element Method 100 where m — 1,2,..., Nst excluding boundary nodes, n = 1,2,..., JVj^ , and i=l,2,... ,M. 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,i0 is the direct capacitance per unit length between conductor i and the reference conductor, and Cmii 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 Ii-tiwC.ndzW- f ) (juC^dzXVi-Vj^Ii + dli (5.10) = l k*i or k = 1 k= 1 FC=1 k?i k^i where Cik is the element in [C] , which is related to the direct capacitances by Cu = c,i0+ f ; cmii (i = i ,2, . . . , io (5.12) fc^i Cik = -Cmih (i,k = l,2,...,K;k^i) (5.13) 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 the field solutions. 5.3.1 The energy method With this method, the electric energy stored in the field under specific boundary con-ditions is calculated, and the elements of [C] are derived from the energy. Once the potential distribution is known, the following equation is used to find the electric energy from the field = \ I cE • Eds = i / eV<p • V<pds i M = (5-15) Chapter 5. Admittance Calculation with Finite Element Method 102 where [^] = (5-16) [<j>Ei] is the node value vector of <f> in element Ei. c^ . is the permittivity in 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 k = i k*i (5.17) 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 2W^ (5.18) By assuming Vi = Vj = V0 ^ 0 and = 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 = \cs,X+1- £ cm,X-\cmtJv02 k=i 1 2 1 K 2 1 2 fc = 1 = \(Cii + Cjj + 2Cij)V02 (5.19) 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, WgF is calculated and Cn is evaluated. By superimposing two solutions, WEp and C t J can then be found. 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 con-ductors 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,<Z2,...,<Zx]r (5-22) 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/dx2+ dy2 (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 Tcr The plus and minus signs in the above formula are related to the direction of D. ID I is calculated from | D | - 4 E | - « | V * | - « ^ 3 = j + ^ (»Y+(»)' (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 = 0, the corresponding integral Aq can be derived as A<? = eEiyJ(xn-xm)* + (yn-ymy £ ±^ (^j + (^j d(n (5.26) and dx dx £.9> \dU dx + SC„ dx. ~du - — d u — P i n { N M ( 5 - 2 9 ) - Pjm(Np,(m) — — (5.30) where Pjm, Pjn, and Pjo are given by (3.5) and (3.6). and in (5-27) and (5.28) can also be expressed in terms of and or in terms of and if (5.29) and (5.30) are modified correspondingly. In order to simplify the calculation of -g^- and Pm(Np,Q is changed into real polynomial form 1 m Pm(Np,0 = —52amiC (5.31) "•m i=0 in which dm is the common denominator. amt- is the polynomial coefficient. dm and amt-are associated with order Np and with m, and they are listed for order 1 to 6 in Tab.5.1. As |£ and | ^ can be found for any given (( m , £„, (0) from (5.27) to (5.30), Aq in (5.26) can now be evaluated by applying the Gaussian quadrature discussed in Section 3.3.1. Chapter 5. Admittance Calculation with Finite Element Method 106 Table 5.1: dm and ami in Pm(Np,Q =.£Efeo«-m<C* m dm Oml <*m2 Ora3 «m4 <&m5 Om6 1 0 1 1 l 1 0 1 2 0 1 1 l 1 0 2 2 1 0 -1 2 3 0 1 1 1 1 0 3 2 2 0 -3 9 3 2 0 2 -9 9 4 0 1 1 1 1 0 4 2 1 0 -2 8 3 3 0 4 -24 32 4 3 0 -3 22 -48 32 5 0 1 1 1 1 0 5 2 2 0 -5 25 3 6 0 10 -75 125 4 24 0 -30 275 -750 625 5 24 0 24 -250 875 -1250 625 6 0 1 1 1 1 0 6 2 1 0 -3 18 3 1 0 2 -18 36 4 2 0 -3 33 -108 108 5 5 0 6 -75 315 -540 324 6 10 0 -10 137 -675 1530 -1620 648 For isoparametric elements, dx, dy, and |£ in element E{ can be derived as dx dy §± r dx —-§± L dy J L dx dx dv dv Sj/ d£ 8y dv - dv (5.32) du du d_ dx (5.33) where [JEi] and [DEi] are given by (3.22) and (3.25), respectively. \0\ is given by (3.19) for 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) The corresponding integral Aq becomes dx dv in dv 1 0 (5.34) (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 dv J J * du L dv i - i (5.36) with A 5 having the same form as in (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. The fcth insulation 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 ek is the permittivity of the fcth insulation. The general form of [C] for SC coaxial cables, which only has three diagonals[10], is c{ cf cf ci Cf 0 where [C} = riod rid riod ° 3 ° 3 ° 4 riod rid riod riod rid cf = C> IN, c? C/Jvt_! + CiNk (k = 2,3,..., K) —CiNu.x (k — 2,3,..., if) (5.38) (5.39) 5.4.2 [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 radii (mm) inner insulation outer insulation division 1 12 18 22 24 division2 12 15 18 22 23 24 division3 12 14 16 18 22 22.7 23.3 24 division4 12 13 15 17 18 22 22.3 23 23.7 24 division5 12 12.3 17.7 18 22 22.1 23.9 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 2 i are the same as — Cu. The capacitances calculated by the surface charge method depend on division scheme. C\2 is the same as — C n and C 2i is slightly different from C12 due to different integration 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 Npth order simplex element, with <j> being a polynomial function of order Np, E calculated from <f> by differentiation will be a polynomial function of order Np — 1. Therefore, the number of sampling points will be ((Np + l)/2), truncated to the nearest integer. For 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 Energy Method (/iF/km) Surface charge method (/iF/km) element Cu (-C121 -C21) Cn (-C12) Cn analytical 0.137017 0.775503 0.137017 -0.137017 0.775503 ISO 0.137037 0.775523 0.133333 -0.133333 0.771014 sim2 0.137039 0.775533 0.133337 -0.133337 0.771034 sim3 0.137018 0.775507 0.136559 -0.137476 0.775724 divisionl sim4 0.137017 0.775505 0.136963 -0.136963 0.775439 sim5 0.137017 0.775504 0.136999 -0.137035 0.775278 sim6 0.137017 0.775504 0.137014 -0.137014 0.775487 ISO 0.137018 0.775504 0.135886 -0.136261 0.774537 siml 0.137503 0.776107 0.122225 -0.150004 0.774416 sim2 0.137019 0.775509 0.135882 -0.136274 0.774382 sim3 0.137017 0.775505 0.136932 -0.137069 0.775334 division2 sim4 0.137017 0.775504 0.137001 -0.137026 0.775290 sim5 0.137017 0.775504 0.137004 -0.137029 0.775312 sim6 0.137017 0.775504 0.137004 -0.137029 0.775346 ISO 0.137017 0.775503 0.136476 -0.136701 0.775082 siml 0.137236 0.775784 0.126680 -0.145309 0.773857 sim2 0.137018 0.775507 0.136467 -0.136719 0.774854 sim3 0.137017 0.775505 0.136978 -0.137045 0.775232 division3 sim4 0.137017 0.775504 0.137000 -0.137035 0.775239 sim5 0.137017 0.775504 0.137000 -0.137036 0.775291 sim6 0.137017 0.775504 0.137000 -0.137036 0.775349 ISO 0.137017 0.775503 0.136871 -0.136942 0.775409 siml 0.137178 0.775717 0.131691 -0.141097 0.775312 sim2 0.137017 0.775505 0.136852 -0.136976 0.774867 division4 sim3 0.137017 0.775504 0.136987 -0.137053 0.774995 sim4 0.137017 0.775504 0.136988 -0.137053 0.775182 sim5 0.137017 0.775504 0.136987 -0.137052 0.775321 sim6 0.137017 0.775504 0.136989 -0.137048 0.775380 ISO 0.137029 0.775515 0.137015 -0.137022 0.775506 siml 0.138374 0.777170 0.136666 -0.139537 0.776884 sim2 0.137029 0.775517 0.136955 -0.137140 0.774046 division5 sim3 0.137017 0.775504 0.136938 -0.137119 0.775028 sim4 0.137017 0.775504 0.136948 -0.137091 0.775356 sim5 0.137017 0.775504 0.136966 -0.137064 0.775324 sim6 0.137017 0.775504 0.136981 -0.137049 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 |E(r)| (V/mm) r=12 mm r=18 mm r=22 mm r=24 mm analytical 0.205525 0.137017 0.522398 0.478865 divisionl 0.200000 0.133333 0.521739 0.478261 division2 0.203829 0.136261 0.522226 0.478720 division3 0.204713 0.136701 0.522312 0.478795 division4 0.205306 0.136943 0.522381 0.478852 division5 0.205522 0.137023 0.522396 0.478863 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 num-bers 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 the figure, 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 112 0.4 s d 0.2 on U 0.1 0.0 EM_iso a = C u (division3) v = C n ( d i v i s i o n s ) -A = C M idiviaion3) x = C M (division5) 20 —r 40 60 0 (degree) 120 0.8 3! o.e-s d 0.4-fe & A 0.2 0.0 -D = C u SCM_iso division31 — O i l (division5) -» = Cn H = Cn A = C M X = C M division3, division5 division3. 60 BO 9 (degree) 120 70 80 H 50 5 40 d • 3 0 u o H 20 w 10 0 EM sim3 - n = C n 7division3) •v = C u (division5) = C M (division3) -x = C u (division5) 20 40 eo eo 8 (degree) 100 120 70 60 50-£i 40 d 30-20 10 H 0 SCM sim3 -a = C u (diviaion3j v = C u (division5) ••o = C»i (division3) -+ = C t i (division5j •A = C M (division3) -x = C M (division5) 120 8 (degree) Figure 5.3: Errors in [C] calculation of a SC coaxial cable for different span angles 5.5 Summary In this chapter general procedures for solving the electrostatic field with FEM are dis-cussed. 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 Sys-tems 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 ca-ble 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 cal-culated 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 air air earth u r, re •s \ s 1 i^=12.1mm '«1=24.5mm h=l.2m j=300mm re=50.9mm fr2=42.3mm rB2=46mm p.= lOOQm LU=un rA3=50.9mm e ^ 0 (a) cable data (b) buried system earth A B C © © © A=1.2m h2 =400mm /j3=100mm S =300mm w =150mm pg=100£2m \le=\iQ (c) tunnel installed 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 [ Z A A ] [ Z A B 1 [ Z A C 1 0]= [ZAB] [ * B B ] [ *BC ] t6-1) . iZAC\ iZBC\ [ZCC\ Subscripts A, B, and C represent the phases. All submatrices are 2x2 matrices, and [•^AA1» [^BB1» a n < ^ [^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. Be-cause the meshes are symmetrical with respect to the center cable, [^AA] =^CCl a n c * [Z^g]=[ZgQ]. With Pollaczek's formula, it is assumed that [^AA]=[^BB] =[^CC3> * n a * all four elements in [^AB ! a r e s a m e > a n a ^ that all f ° u r elements in [^ACl a r e ^ n e s a m e -Tab.6.1 and Tab.6.2 confirm that these assumptions are reasonably accurate for the given pe=100flm. The tables show that the differences between the results from Pollaczek's 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 117 Table 6.1: [ • Z A A K I ^ C C ' D [ ^ B B I °^ *^ e buried three-phase cable system / (Hz) R (ft/km) L (mH/km) Rn R12 R22 £11 L12 L22 60 [*AA1 conv. prop. 0.0748300 0.0597006 0.262705 0.0749538 0.0598244 0.262829 2.10829 1.96073 1.95793 2.11098 1.96342 1.96062 [*BB! conv. prop. 0.0750753 0.0599458 0.262960 0.0751972 0.0600678 0.263072 2.10904 1.96147 1.95868 2.11172 1.96416 1.96137 Pollaczek 0.0745183 0.0594093 0.262382 2.12015 1.97267 1.96988 6000 [*AA! conv. prop. 6.37885 6.16838 6.33678 6.39170 6.18123 6.34963 1.61213 1.49146 1.48897 1.61402 1.49335 1.49086 [*BB1 conv. prop. 6.38720 6.17672 6.34512 6.39999 6.18952 6.35792 1.61034 1.48967 1.48718 1.61225 1.49158 1.48909 Pollaczek 6.37657 6.16547 6.33386 1.62790 1.50717 1.50468 600000 [*AA! conv. prop. 694.949 691.010 691.011 696.979 693.041 693.041 1.09676 0.98654 0.98654 1.09769 0.98747 0.98747 t*BBl conv. prop. 694.819 690.880 690.881 696.919 692.980 692.981 1.09482 0.98460 0.98460 1.09576 0.98554 0.98554 Pollaczek 698.287 694.327 694.327 1.11182 1.00154 1.00154 Table 6.2: [ ^ A B K I ^ B C J ) t^ACl °^ buried three-phase cable system / (Hz) J2 (ft/km) L (mH/km) Ru R12 -R21 R22 i l l Ll2 L21 L22 60 [*AB1 conv. prop. Poll. 0.05944 0.05944 0.05944 0.05944 0.05957 0.05957 0.05957 0.05957 0.05940 1.586 1.586 1.586 1.586 1.589 1.589 1.589 1.589 1.589 f*ACl conv. prop. Poll. 0.05895 0.05895 0.05895 0.05895 0.05908 0.05908 0.05908 0.05908 0.05940 1.449 1.449 1.449 1.449 1.452 1.452 1.452 1.452 1.451 6000 [*AB1 conv. prop. Poll. 6.085 6.085 6.085 6.085 6.098 6.098 6.098 6.098 6.092 1.120 1.120 1.120 1.120 1.122 1.122 1.122 1.122 1.125 f*ACl conv. prop. Poll. 6.066 6.066 6.066 6.066 6.079 6.079 6.079 6.079 6.091 0.989 0.989 0.989 0.989 0.990 0.990 0.990 0.990 0.986 600000 P A B I conv. prop. Poll. 685.9 685.9 685.9 685.9 687.9 687.9 687.9 687.9 689.1 0.622 0.622 0.622 0.622 0.623 0.623 0.623 0.623 0.626 [*ACl conv. prop. Poll. 679.9 679.9 679.9 679.9 681.8 681.8 681.8 681.8 682.5 0.492 0.492 0.492 0.492 0.492 0.492 0.492 0.492 0.488 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. Case Studies in [Z] and [C] Calculations 118 Tab.6.4. re=300mm is used with Pollaczek's formula. Again, the differences between Pollaczek's formula and the conventional FEM become noticeable at high frequencies. Table 6.3: [^AAK^CC]) a n < * I ^ B B I °^ ^ e * u n n e l installed three-phase cable system / (Hz) # (ft/km) L (mH/km) Rll #12 #22 L\\ Lu L22 60 1*AA1 conv. prop. 0.0748252 0.0596958 0.262700 0.0749517 0.0598223 0.262827 2.11131 1.96375 1.96096 2.11394 1.96637 1.96358 [*BB1 conv. prop. 0.0750711 0.0599417 0.262946 0.0751999 0.0600705 0.263075 2.11116 1.96360 1.96080 2.11376 1.96620 1.96340 Pollaczek 0.0745182 0.0594092 0.262382 2.12015 1.97267 1.96988 6000 1*AA1 conv. prop. 6.37058 6.16010 6.32850 6.38852 6.17804 6.34644 1.61540 1.49473 1.49224 1.61716 1.49649 1.49400 [*BB1 conv. prop. 6.37804 6.16757 6.33597 6.39624 6.18576 6.35416 1.61276 1.49209 1.48360 1.61448 1.49381 1.49132 Pollaczek 6.37618 6.16508 6.33347 1.62790 1.50717 1.50468 600000 [*AA1 conv. prop. 675.093 671.155 671.155 685.630 681.692 681.692 1.10341 0.99319 0.99319 1.10176 0.99154 0.99154 [*BB! conv. prop. 672.579 668.641 668.641 683.822 679.884 679.884 1.10084 0.99061 0.99061 1.09903 0.98881 0.98881 Pollaczek 695.484 691.523 691.523 1.11193 1.00166 1.00166 Table 6.4: [ - ^ A B K I ^ B C ^ ) a n c ^ l^ACl °^ tunnel installed three-phase cable system / (Hz) # (ft/km) L (mH/km) #11 #12 #21 #22 Lu L12 L21 L22 6 0 P A B I conv. prop. Poll. 0.05944 0.05944 0.05944 0.05944 0.05957 0.05957 0.05957 0.05957 0.05940 1.584 1.584 1.584 1.584 1.587 1.587 1.587 1.587 1.589 [*ACl conv. prop. Poll. 0.05894 0.05894 0.05894 0.05804 0.05907 0.05907 0.05907 0.05907 0.05940 1.448 1.448 1.448 1.448 1.450 1.450 1.450 1.450 1.451 6 0 0 0 I*AB1 conv. prop. Poll. 6.077 6.077 6.077 6.077 6.095 6.095 6.095 6.095 6.092 1.118 1.118 1.118 1.118 1.120 1.120 1.120 1.120 1.125 [*ACl conv. prop. Poll. 6.060 6.060 6.060 6.060 6.078 6.078 6.078 6.078 6.091 0.987 0.987 0.087 0.987 0.989 0.989 0.989 0.989 0.986 6 0 0 0 0 0 [*AB! conv. prop. Poll. 665.5 665.5 665.5 665.5 676.0 676.0 676.0 676.0 689.1 0.624 0.624 0.624 0.624 0.622 0.622 0.622 0.622 0.626 t^ ACl conv. prop. Poll. 661.4 661.4 661.4 661.4 671.9 671.9 671.9 671.9 682.5 0.494 0.494 0.494 0.494 0.492 0.492 0.492 0.492 0.488 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 PT Cables 6.3.1 The [Z] calculation of PT 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 JRi = 6903.9 - 47.4-0'-1) hi = 6928.1 - 47.2 < i -1) /«, =6760.5 - 44.3 •(/-!) A* = 6857.8 - 41.7.(1-, =36453-1340«0-1) =36486 - 1339'(j-1) //, =31626-1105-(. -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 the field solutions were derived for a loop current between a single current filament and the pipe, the influence of the other non-current carrying conductors on the field distributions inside the pipe and within the pipe conduc-tor is ignored. The influence of the finite dimension of the current carrying conductor, which is replaced by a current filament for 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 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 240 mm (a) cable data (b) triangle arrangement (c) cradle arrangement Figure 6.4: A 230kV PT cable system 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 con-ductor. 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 magnetic flux lines[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 123 Table 6.5: [ - ^ A I A N < ^ [ ^ B B K ^ C C " ] ) °^ c a D ^ e with a triangle arrangement / (Hz) # (n/km) L (mH/km) #11 #12 #22 I'll L\2 2/22 60 I*AA1 FEM ana. 0.282226 0.260819 0.662920 0.275526 0.254153 0.656159 1.01507 0.86707 0.86545 1.03926 0.89132 0.88969 [%B! FEM ana. 0.285957 0.264551 0.666666 0.279190 0.257818 0.659823 1.04822 0.90023 0.89862 1.07678 0.92884 0.92721 600 1*AA1 FEM ana. 0.975600 0.916694 1.31873 0.908456 0.849402 1.25125 0.48853 0.36789 0.36628 0.53731 0.41660 0.41498 [ZBBI FEM ana. 1.01266 0.053773 1.35584 0.977012 0.917958 1.31981 0.50442 0.38378 0.38216 0.55370 0.43298 0.43137 6000 P A A I FEM ana. 3.20001 3.00974 3.40539 3.20516 3.01413 3.40962 0.29346 0.18268 0.18108 0.36005 0.24922 0.24762 13BB] FEM ana. 3.53022 3.33996 3.73563 3.89726 3.70623 4.10172 0.29957 0.18879 0.18719 0.34546 0.23464 0.23304 60000 [*AA! FEM ana. 11.9585 10.6357 10.7013 12.4560 11.1361 11.2030 0.23325 0.12740 0.12666 0.29353 0.18764 0.18691 (ZBB1 FEM ana. 16.3469 15.0241 15.0897 16.3205 15.0007 15.0675 0.22250 0.11665 0.11591 0.25717 0.15128 0.15055 600000 I*AAJ FEM ana. 42.6742 38.1846 38.1927 41.3592 36.8655 36.8741 0.20765 0.10516 0.10516 0.26805 0.16552 0.16552 I^BB) FEM ana. 76.1295 71.6390 71.6480 56.6433 52.1496 52.1582 0.18052 0.07803 0.07803 0.22257 0.12004 0.12004 Table 6.6: [ ^ A B K I ^ A C " ] ) [^ BCl °^ ^ e PT cable with a triangle arrangement / (Hz) # (ft/km) L (mH/km) #11 #12 #21 #22 Lu 2/12 2/2i 2)2 60 I*AB1 FEM ana. 0.2498 0.2498 0.2498 0.2498 0.2493 0.2493 0.2493 0.2493 0.6857 0.6857 0.6857 0.6858 0.6824 0.6824 0.6824 0.6824 [*BCl FEM ana. 0.2470 0.2470 0.2470 0.2470 0.2482 0.2482 0.2482 0.2482 0.6470 0.6470 0.6470 0.6470 0.6396 0.6396 0.6396 0.6396 600 I*AB1 FEM ana. 0.7763 0.7764 0.7763 0.7764 0.7611 0.7611 0.7611 0.7611 0.2357 0.2357 0.2357 0.2357 0.2357 0.2357 0.2357 0.2357 [*BC] FEM ana. 0.7529 0.7528 0.7528 0.7528 0.7425 0.7425 0.7425 0.7425 0.2055 0.2055 0.2055 0.2055 0.1988 0.1988 0.1988 0.1988 6000 1*AB1 FEM ana. 2.506 2.506 2.506 2.506 2.184 2.184 2.184 2.184 0.0947 0.0947 0.0947 0.0947 0.1069 0.1069 0.1069 0.1069 [*BCl FEM ana. 2.372 2.372 2.372 2.372 2.010 2.010 2.010 2.010 0.0722 0.0722 0.0722 0.0722 0.0782 0.0782 0.0782 0.0782 60000 [*AB! FEM ana. 8.729 8.729 8.729 8.729 6.283 6.283 6.283 6.283 0.0466 0.0466 0.0466 0.0466 0.0711 0.0711 0.0711 0.0711 [*BCl FEM ana. 7.696 7.696 7.696 7.696 5.428 5.428 5.428 5.428 0.0282 0.0282 0.0282 0.0282 0.0474 0.0474 0.0474 0.0474 600000 1*AB1 FEM ana. 33.58 33.58 33.58 33.58 18.98 18.98 18.98 18.98 0.0283 0.0283 0.0283 0.0283 0.0604 0.0604 0.0604 0.0604 [*BCl FEM ana. 25.02 25.02 25.02 25.02 15.86 15.86 15.86 15.86 0.0136 0.0136 0.0136 0.0136 0.0386 0.0386 0.0386 0.0386 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 I 16.0 H 15.0 14.04 13.0 12.0 -° = one cable v = three cables 0.030 —i i 1 1 1 1 1— 4S 80 135 180 225 270 315 380 0 (degrees) 0.025 & 0.020 o.oisH 0.010H 0.005 -o = one cable — i 1 1 1 1 1 1— 0 45 80 135 180 225 270 315 380 <t> (degrees) (a) surface \ J\ at 60Hz (b) surface \ J\ at 6kHz Figure 6.7: \J\ on the inner surface of the pipe caused by the loop current 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, rJi2=20mm, rBj=23mm, and rJ43=25mm for the PT cable in Fig. 6.4(b) with 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 PT 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 126 Table 6.7: [ ^ A A K ^ C C D [ ^ B B I °f the PT cable with a cradle arrangement / (Hz) # (fi/km) L (mH/km) #11 #12 #22 Xll X12 L22 60 1*AA1 FEM ana. 0.285638 0.264232 0.666324 0.279190 0.257818 0.659823 1.04983 0.90185 0.90024 1.07678 0.92884 0.92721 I^BBI PEM ana. 0.287547 0.266144 0.668125 0.279190 0.257818 0.659823 1.04277 0.80478 0.80317 1.07678 0.92884 0.92721 600 P A A I FEM ana. 1.01671 0.957819 1.35986 0.977012 0.917958 1.31981 0.50716 0.38652 0.38491 0.55370 0.43298 0.43137 [^BB! PEM ana. 1.04002 0.981136 1.38306 0.977012 0.917958 1.31981 0.40151 0.37087 0.36026 0.55370 0.43298 0.43137 6000 I*AA1 FEM ana. 3.65653 3.46628 3.86192 3.89726 3.70623 4.10172 0.29510 0.18432 0.18272 0.34546 0.23464 0.23304 (^BB) FEM ana. 3.56767 3.37744 3.77298 3.89726 3.70623 4.10172 0.27689 0.16611 0.16451 0.34546 0.23464 0.23304 60000 [*AA1 FEM ana. 16.2286 14.9059 14.9714 16.3205 15.0007 15.0675 0.21712 0.11126 0.11053 0.25717 0.15128 0.15055 I^BB) FEM ana. 14.9878 13.6651 13.7305 16.3205 15.0007 15.0675 0.20428 0.09843 0.09769 0.25717 0.15128 0.15055 600000 t^ AAl FEM ana. 73.0357 68.5463 68.5544 56.6433 52.1496 52.1582 0.17656 0.07407 0.07407 0.22257 0.12004 0.12004 [%B1 FEM ana. 64.4238 50.9340 50.0430 56.6433 52.1496 52.1582 0.16810 0.06570 0.06570 0.22257 0.12004 0.12004 Table 6.8: [ ^ A B K ^ B C ] ) a n c ^ l^ACl °^ PT cable with a cradle arrangement / (Hz) # (ft/km) L (mH/km) #11 #12 #21 #22 £ n L\2 L21 L22 60 P A B I FEM ana. 0.2513 0.2513 0.2513 0.2513 0.2508 0.2508 0.2508 0.2508 0.7003 0.7003 0.7003 0.7003 0.6981 0.6981 0.6981 0.6981 l*ACl PEM ana. 0.2446 0.2446 0.2446 0.2446 0.2461 0.2461 0.2461 0.2461 0.6115 0.6115 0.6115 0.6115 0.5975 0.5975 0.5975 0.5975 600 1*AB1 FEM ana. 0.7972 0.7972 0.7972 0.7972 0.7855 0.7855 0.7855 0.7855 0.2426 0.2426 0.2426 0.2426 0.2429 0.2429 0.2429 0.2429 [*ACl FEM ana. 0.7149 0.7149 0.7149 0.7149 0.7071 0.7071 0.7071 0.7071 0.1845 0.1845 0.1845 0.1845 0.1689 0.1689 0.1689 0.1689 6000 [*AB! FEM ana. 2.663 2.663 2.663 2.663 2.351 2.351 2.351 2.351 0.0940 0.0940 0.0940 0.0940 0.1046 0.1046 0.1046 0.1046 [*ACl FEM ana. 2.218 2.218 2.218 2.218 1.759 1.759 1.759 1.759 0.0635 0.0635 0.0635 0.0635 0.0622 0.0622 0.0622 0.0622 60000 [*AB! FEM ana. 9.605 9.605 9.605 9.605 6.832 6.832 6.832 6.832 0.0407 0.0407 0.0407 0.0407 0.0651 0.0651 0.0651 0.0651 [*AC] FEM ana. 7.024 7.024 7.024 7.024 4.502 4.502 4.502 4.502 0.0218 0.0218 0.0218 0.0218 0.0374 0.0374 0.0374 0.0374 600000 1*AB1 FEM ana. 35.39 35.39 35.39 35.39 20.51 20.51 20.51 20.51 0.0208 0.0208 0.0208 0.0208 0.0534 0.0534 0.0534 0.0534 t*ACl FEM ana. 21.59 21.50 21.50 21.59 12.93 12.93 12.93 12.93 0.0094 0.0094 0.0094 0.0094 0.0302 0.0302 0.0302 0.0302 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.243091 0 -0.040218 0 -0.040218 0 0 0.109652 -0.109652 0 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.109442 0 0 0 0 -0.109592 0.243030 0 -0.040217 0 -0.040217 0 0 0.109442 -0.109442 0 0 0 -0.040217 -0.109593 0.341983 0 -0.014692 0 0 0 0 0.109442 -0.109442 0 -0.040217 0 -0.014692 -0.109593 0.341988 [C] found with an alytical formulas 0.109795 -0.109795 0 0 0 0 -0.109795 0.209431 0 -0.039429 0 -0.039429 0 0 0.109795 -0.109795 0 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.348885 0 -0.053153 0 -0.007603 0 0 0.109651 -0.109651 0 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.109442 0 0 0 0 -0.109593 0.348820 0 -0.053152 0 -0.007603 0 0 0.109440 -0.109440 0 0 0 -0.053152 -0.109590 0.386443 0 -0.053165 0 0 0 0 0.109442 -0.109442 0 -0.007603 0 -0.053165 -0.109593 0.348831 [C] found with an alytical formulas 0.109795 -0.109795 0 0 0 0 -0.109795 0.243352 0 -0.057283 0 -0.007962 0 0 0.109795 -0.109795 0 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 128 6.4 [Z] and [C] Calculations of Sector-Shaped Cables 6.4.1 The [Z] calculation of sector-shaped cables 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 rB = TJ— (6-2) 27T inner radius = )jrB — (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=300mm2 and Lc=70.834mm, with the equivalent radii be-ing rj4=5.622mm and rB=11.274mm. It is assumed that the three equivalent circular 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 r4 in Fig. 6.8(a). The FEM mesh around the cable at 60kHz is plotted in Fig. 6.8(b). The impedances of the cable are listed in Tab.6.11. Due to symmetry, Zi 1=^22=^33 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 / (Hz) R (Jtykm) L (/iH/km) Rn R12 6 FEM app 2.84006 2.78246 2.83996 2.78239 232.020 40.4454 196.081 25.0563 60 FEM app 2.84987 2.78317 2.84334 2.78105 220.673 40.1914 191.259 27.3958 600 FEM app 2.96756 2.78308 2.93723 2.75504 156.797 35.7225 164.736 36.4771 6000 FEM app 3.52505 2.69462 3.52328 2.58640 120.400 38.3245 125.467 47.2753 60000 FEM app 5.15051 2.88824 5.89100 2.22645 105.090 40.6921 108.540 52.1332 600000 FEM app 15.9166 8.82314 18.8742 6.42831 99.0031 38.0609 100.191 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 magnetic field is 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 Ln 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 er=l. If the potential of the upper conductor is IV 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 ^33) Cl2(Cl3>(723) FEM (energy method) 0.147049 -0.0401938 FEM (surface charge method) 0.147323 -0.0402814 analytical formula 0.181528 -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 re-sistance 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 rea = r.(0.92122679 - 1 (6.5) e* e V 4.3856939n2 - 2.3071869n -1.2479854^ K } where r e is the outer radius of the stranded conductor, n is the same as in (6.4). This 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 cal-culation 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. Case Studies in [Z] and [C] Calculations 133 real part of the impedance for such a single conductor system would be the internal re-sistance 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 max-imum 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) Rc (A/km) FEM GSW formula subdivision method Borges da Silva's formula #,=2.25 #,=1.6 43 1.0446 1.4524 1.0328 1.3479 1.0215 80 1.4160 1.9811 1.4088 1.6305 1.3864 100 1.5905 2.2149 1.5751 1.8187 1.5479 130 1.7951 2.5254 1.7959 2.1126 1.7622 max diff with FEM 40.68% 1.13% 29.04% 2.68% given at the bottom of Tab.6.13. With re=13mm and n=12, re,=11.7171mm is found 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) \Jt\ = 3520.7- (i-l) (a) |J | at 60Hz (A/m2) (b) \J\ at 60kHz (A/m2) 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, re=7.8mm, and re,=6.8578mm. The differences Table 6.14: The internal resistance of the one-layer stranded conductor / (kHz) Rc (n/km) FEM GSW formula subdivision method Borges da Silva's formula Kf=2.25 #,=1.6 20 1.2436 1.7335 1.2327 1.7122 1.2246 43 1.7882 2.5418 1.8075 1.8466 1.7688 80 2.4094 3.4669 2.4654 2.4815 2.3920 100 2.7013 3.8761 2.7564 2.8160 2.6677 130 3.0466 4.4195 3.1427 3.2965 3.0338 max diff with FEM 45.06% 3.15% 37.68% 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 re will be 18.2mm and 23.4mm, respectively, with req= 16.5286mm Chapter 6. Case Studies in [Z] and [C] Calculations 135 for the three-layer stranded conductor and req= 21.3293mm for the four-layer stranded 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) Rc (O/km) FEM GSW formula Borges da Silva's formula #,=2.25 #,=1.6 20 0.50658 0.69338 0.49307 0.49425 60 0.87271 1.2010 0.85403 0.84906 100 1.1174 1.5505 1.1025 1.0934 140 1.3169 1.8345 1.3045 1.2920 180 1.4950 2.0802 1.4792 1.4637 220 1.6610 2.2997 1.6353 1.6172 max diff with FEM 39.30% 2.67% 2.71% Table 6.16: The internal resistance of the four-layer stranded conductor / (kHz) Rc (n/km) FEM GSW formula Borges da Silva's formula tf/=2.25 If) =1.6 20 0.39099 0.53337 0.37929 0.38133 60 0.67642 0.92383 0.65694 0.65630 100 0.86673 1.1927 0.84811 0.84563 140 1.0229 1.4112 1.0035 0.99951 180 1.1627 1.6001 1.1379 1.1326 220 1.2934 1.7690 1.2580 1.2515 max diff with FEM 37.95% 2.99% 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 to field tests[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 reason-ably accurate self and mutual impedances for multiphase buried cables, and for tunnel installed cables using approximate inner earth radii. Good agreement between the pro-posed 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 conduc-tor 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 fre-quency 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 resis-tance 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 agree-ment 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, the finite element 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 electric fields were 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] cal-culation from the field solution 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 complete field solution is not required. 2. Quadratic isoparametric elements proved to be more efficient than high-order sim-plex 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 con-ventional FEM if the field truncation boundary is at least ZSe away from the cables {n > 36e), and if the earth is divided with the pattern 10n, 10n+s, 10N +3, 10n + 1 in each decade in the radial direction. For shallowly buried SC cables, the iterative results showed that rj > 12Se is required to make the differences in [Z] from two 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 />e=100f2m and with re varying between 24mm and lm, the maximum errors with 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 shal-lowly buried SC cable in Chapter 4 showed that accurate results were obtained with Pollaczek's formula when re/Se is small. For re=24mm, the maximum differ-ences between the results from Pollaczek's formula and from the FEM are less than 1% in [R] if re/6e <0.03 and in [L] if re/8e <0.095. For large re/«5e, the differences become large. With typical ranges of pe and r e, however, the maximum differences between the two approaches from 1Hz to 1MHz are reasonably small. With pe varying between lOOOfim to l£2m and with r e varying between 24mm to 250mm, 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 e is used. For typical earth resis-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 magnetic flux distribution in the pipe and the current density-distribution 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 frequen-cies. 7. For a sector-shaped cable with a magnetic sheath, the the analytical formulas sug-gested by Ametani produce reasonably accurate resistances and inductances in the low frequency range. If the sheath is non-magnetic or nonexistent, then the in-ductance 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 sim-ple 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 con-sidered 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 ap-proximate 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 Determin-ing 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 Power-Transmission 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 Ef-fect 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 Cal-culating 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 Infinitely-Long Single-Conductor Lines by the Complete Field Solution Method," Proc. IEE, References 143 Vol.125, June 1978, pp.511-517. [20] S. Cristina and M. D'Amore, "Propagation on Polyphase Lossy Lines: A New Pa-rameter Sensitivity Model," IEEE Trans. PAS, vol.PAS-98, Jan./Feb. 1979, pp.35-44. [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 Induc-tive Skin Effect in Rectangular Conductors," IBM J. RES. DEVELOP, vol.23, Nov. 1979, pP652-660. [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 Electro-magnetic 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 Mul-ticonductor 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 Electro-magnetics," IEEE Trans. MAG, vol.MAG-21, Sep. 1985, pp.1829-1834. [36] J. Poltz and E. Kuffel, "Application of Boundary Element Techniques for 2D Eddy-Current 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 Frequency-Dependent 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 Investiga-tion 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 Ca-ble 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 Mag-netics, 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. 2 symm. 0 1 1 2 0 - 1 1 1 1 2 Table A.2: [QW] and [Ts] of the 2nd order simplex element (a) (b) [TS] common denominator = 6 common denominator = 180 0 6 0 8 symm. 0 32 symm. 0 -8 8 0 16 32 0 0 0 3 -1 0 -4 6 0 0 0 -4 8 -4 16 16 0 32 0 0 0 1 -4 3 -1 -4 0 -1 0 6 146 Appendix A. Integral matrices [QW] and [Ts] of simplex elements Table A.3: [QM] and [Ts] of the 3rd order simplex element (a) [Q(1)] common denominator = 80 (b) [Ts] common denominator = 6720 0 76 0 135 symm. 18 540 symm. 0 -135 135 18 270 540 0 -27 27 135 0 -189 -135 540 0 0 0 -162 324 36 162 162 162 1944 0 27 -27 27 -162 135 0 -135 -189 -54 162 540 0 3 -3 3 0 -3 34 11 0 27 18 36 27 76 0 0 0 0 0 0 -54 135 27 -135 -54 270 162 -135 18 540 0 0 0 0 0 0 27 -108 135 27 -54 -135 -135 162 270 0 -189 0 -3 3 -3 0 3 -7 27 -54 34 11 27 0 27 36 18 11 0 Table A.4: [QW] and [Ts] of the 4th order simplex element (a) [<2(1)] (common denominator = 1890) 0 0 3968 0 -3968 3968 0 -1440 1440 4632 0 0 0 -5376 10752 symm, 0 1440 -1440 744 -5376 4632 0 640 -640 -1248 1536 -288 3456 0 0 0 768 -1536 768 -4608 10752 0 0 0 768 -1536 768 1536 -7680 10752 0 -640 640 -288 1536 -1248 -384 1536 -4608 3456 0 -80 80 80 -160 80 240 -160 -160 80 705 0 0 0 -128 256 -128 -256 256 256 -256 -1232 3456 0 0 0 96 -192 96 192 -192 -192 192 884 -3680 0 0 0 -128 256 -128 -256 256 256 -256 -464 1920 0 80 -80 80 -160 80 80 -160 -160 240 107 -464 5592 -3680 3456 884 -1232 705 (b) [Ts] (common denominator = 56700) 290 160 2560 160 1280 2560 -80 -1280 -960 3168 symm. 160 1280 1280 384 10752 -80 -960 -1280 48 384 3168 0 768 512 -1280 256 64 2560 -160 256 -256 384 -1536 -768 1280 10752 -160 -256 256 -768 -1536 384 -256 -1536 10752 0 512 768 64 256 -1280 256 -256 1280 2560 -27 0 -112 -80 -160 -12 160 160 -160 -112 290 -112 512 256 -960 -256 64 1280 1280 256 512 160 2560 -12 64 64 48 -768 48 -960 384 384 -960 -80 -1280 3168 -112 256 512 64 -256 -960 512 256 1280 1280 0 768 -1280 -27 -112 0 -12 -160 -80 -112 -160 160 160 -27 0 -80 2560 160 290 (a) (common denominator = 290304) 0 743750 0 -743750 743750 0 -405000 405000 1072500 0 0 0-1150000 2300000 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 0 -150000 300000 -150000 -100000 412500 -525000 212500 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 0 -5000 10000 -5000 30000 -92500 95000 -32500 0 0 0 -5000 10000 -5000 -325O0 95000 -92500 30000 0 0 0 20625 -41250 20625 32500 -79375 61250 -14375 0 -11850 11850 -11025 31250 -20225 -11025 31250 -15625 -4600 (b) [771 (commoi 53244 45270 563500 45270 281750 563500 -36720 -367250 -281750 846000 34200 392500 392500 -55000 2700000 -36720 -281750 -367250 113000 -55000 846000 17730 295250 220750 -516500 195000 -52000 846000 -37350 3750 -146250 202500 -600000 -210000 202500 2925000 -37350 -146250 3750 -210000 -600000 202500 -22500 -450000 2925000 17730 220750 295250 -52000 195000 -516500 108000 -22500 202500 846000 -880 -152125 -95125 295250 -107500 -1750 -367250 3750 -52500 -55500 33700 -107500 30000 195000 450000 120000 -55000 -600000 150000 120000 24900 -52500 -52500 -22500 150000 -22500 •210000 -450000 -450000 -210000 33700 30000 -107500 120000 450000 195000 120000 150000 -600000 -55000 -880 -95125 -152125 -1750 -107500 295250 -55500 -52500 3750 •367250 4747 -880 13935 17730 33700 7940 -36720 -37350 24900 7940 13935 -95125 -15750 220750 30000 -55500 -281750 -146250 -52500 -1750 7940 -1750 -55500 -52000 120000 108000 113000 -210000 -22500 -52000 7940 -55500 -1750 108000 120000 -52000 -52000 -22500 -210000 113000 13935 -15750 -95125 -55500 30000 220750 -1750 -52500 -146250 -281750 4747 13935 -880 7940 33700 17730 7940 24900 -37350 -36720 symm. 577500 -818750 2050000 393750-1837500 2887500 -193750 800000-1837500 2050000 41250 -193750 393750 -818750 577500 58725 -62500 -15625 31250 -11850 99402 -83125 146250 14375 -135000 57500 -188125 577500 83125 -177500 1250 197500 -104375 172500 -745625 1282500 -104375 197500 1250 -177500 83125 -130625 570000-1148750 1282500 575O0 -135000 14375 146250 -83125 58750 -272500 570000 -745625 577500 -11850 31250 -15625 -62500 58725 -11902 58750 -130625 172500 -188125 99402 19160064) symm. 392500 2700000 -146250 -600000 2925000 30000 450000 -600000 2700000 -15750 30000 -146250 392500 563500 45270 34200 -37350 33700 13935 53244 281750 392500 3750 -107500 -95125 45270 -281750 -55000 202500 195000 220750 -36720 220750 195000 202500 -55000 -281750 17730 -95125 -107500 3750 392500 281750 -880 13935 33700 -37350 34200 45270 4747 563500 -367250 846000 295250 -516500 •152125 295250 -880 17730 0 I—* n > Cn 55 fcr* a Cn XT O H &* n B *—• n X £L n B n P 846000 •367250 563500 -36720 45270 53244 Appendix A. Integral matrices [Q^] and jTs] of simplex elements Table A.6: [QW] and [Ts] of the 6th order simplex element 149 8 ON I-i O «-» c3 C s o c <o T3 C o o o & E S S z s M i l l : 2 S G ' *9 T ^ 9 i l l i l l » ? n lis IH s a a - 8 : go r> lis S J 3 H S SS 3 2 8 5 3 J S S i Of W N « o o r« so e S 8 X 8 3 S 8 a § 2 S 3 3 S 5 ^ « A i cn* *f 3 I s i s ' ? ! 9 J S 8 I s a 1 n o o vo S S 88 S ' saps rt M j j i e rt oe s ? s 5! R M a s 3 f 8 8 3 r — cfi «n S3 S tsSI s s a « N 09 V « N N R 5 3 $ a s a " 2 S S S S 2 8 3 3 3 s s a s R I R K « a R !S * <* - s 5 : I 5 8 s a * s I 8 2 3 S S 8 8 S S S ! i § s p; r« <n • s s a s lap S § 1 s s a s s a « s sp a a 1 s ! : s i 33 P 8 a s; a § s : o3 S S 8 a s i '•as I S O 53 8 X s p a list §3 a s s 223 te ^ a • 2 a T s s p s s ^ s a ^ a a c ? a s s x m tn oo & SSSR-T T s s <»"i ft yj «o 2 a a § s | l a s S P £ !8 I m o « «o R S1S I § s i 2 a: o so r* s a R a a 2 a a S R •* £ J£ IIS s ° ° ° a R £5 a a: © O O ' » * 8 0 0 8 S 0 O O a R o o o o o o o o o o o o o o o e o o o o s a' • ° ° a R si a O O O O I o s e o c •3 c o o S § s III I 8 8 8 i l l ? l i ii SSI3 9 « * i 5 i f f l i s 3 a till Mil S S I g i I i 1111 8 a g « : a S I S : s a B P i s g I 8 | x = g s 8 I I g 9 « o » « 1111 s M * 8 a x a i l l s - j s ? s a a j si?** 1 5 1 S I l l s s s | a H I 8 B 5 i l = 5 tt « 5 8 a g 5 g T g i i £§ I ' M 1111! 8 S 8 f i i 11 V R a a § § 1 J I I c 11 11 ' - ^ T « i 1 p s s * s T S 3 T = S S s a p P C S * : 5 s -8 8 § i m m | i s 11 § 8 I S 8 § 8 i l l a S 5 a a a 111 | 8 11 " l i a 5 "? S s »* 11 8 a s a S3 15 S 8 0 X l i s — r> -» § 1 1 C S S I 1 l i n ? i § I 5§ M 3 8 O O O O 1 i o o o o s a a a a * » R • 2 i i 3 2 S i ?5 E 11 S g 2 1 3 a B 5 s § s a 11 s i $ ? s"? P S ! § s I § 8 51 S "7 !§ ; R l i s a s a P SiSe 11 If: 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. O region 1: air (M-fl) (xryf) sf current filament • ,( x f,tf)m Figure B.l: A current filament buried in the earth -I \ 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 x B = J V E = 0 (B.2) (B.3) 150 Appendix B. Detailed Derivation of Pollaczek's Formula 151 the following equation is got V x V x E = V(V • E) — V 2 E = - V 2 E = -jwV x B = -jwpJ = ^ ^ E (B.4) P The current density has one direction only. Therefore, the two-dimensional principal equation descibing the field in the problem is V2E = jvfiJ = \ E (B.5) P where P2 = -A- (B.6) p is the complex penetration depth. If EA and EE respectively represent E fields in the air and in the earth, (B.5) can now be splitted into two equations V2EA = 0 y>0 (B.7) V2EE = \E. + jup.I6(x-xf)6(t-t,) t>0 (B.8) Pi w here Pe = "T~~ (B.9) jup,e I is the magnitude of the filament current. The Delta function is associated with the filament. The boundary conditions are Ea = Ee = E0 y = 0 (B.10) ldEa ldEe ldEe = — = — y = 0; t = 0 (B.ll) p,a ay p,e ay p,e at dEa dEe dEa dEe n a = e = ~^ = ~dT = ~dV = ^ ' = * = 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) = / f{x)e~3U*dx (B.13) J—oo Appendix B. Detailed Derivation of Pollaczek's Formula 152 /(*) = ^ J~J(a)e>xada (B.14) To y and t the following Fourier sine transformation pair is applied 71(A) = [°° f(y)sm(Xy)dy Jo f(y) = - rjt{\)sm{y\)d\ 7T JO For region 1, x and y are used as the coordinates. Define K = [°° Eae-ja*dx J—oo Ea = / Easin(Xy)dy Jo Applying the Fourier transformation to x as r°° d2E I e~3axdx J-oo dx2 = ( - J « ) dEa dx -jax /°° 8E oo oo = - Eae-iax(-ja)2 °° + f°° Eae-jax(-ja)2dx -oo J-00 /oo Eae~}axdx = -a2Ea -00 r d2Ea By2 e~3axdx = Then (B.7) becomes d2E~a dy2 Applying the Fourier sine transformation to y as Jroo ' -a2 Eas\n{Xy)dy = -a2Ea 0 d2E2 . dWa . svn\<Xy)dy = -—— sin(Ay) 1° Jo dy2 dy f°° dE dy .oo roo = - AjEocos(Ay)|0 - A 2 y Easm(\y)dy = XEQ — X2EA (B.15) (B.16) (B.17) (B.18) (B.19) (B.20) (B.21) (B.22) (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 r<x> Ee = / Eee~3axdx (B.25) J — OO We = rE~esm(\t)dt (B.26) Jo By similar process, the following equations can be derived -a2Te + ^ = ±--Ee+jupJe-ia*'8{t-tf) (B.27) dy2 pi {a2 + \ + \2)Te = \E0'-jupeIe-jax'Sm{\tf) (B.28) pi By assigning 02 = a2 + \ (B.29) Pi The inverse Fourier sine transformation is applied to (B.24) and (B.30) to get Ea and K = I* (mi) According to the mathematical handbook f - ^ m d X . * - a > „ ; y > 0 (B.33) Jo a* + M 2 y>0 (B.34) (B.28) becomes 7T 2' r S i n ( y ) S ! f A)<*A = JL («- (« / - ' ) •_ *>0; *,>*><> (B.35) Jo 02 + A2 40 v ' = i L (e-(t-t/)# _ e-(*+t/)«) 0>fj; *>*,>() (B.36) = J L ( e-|t-t/ll«l _ e-lW/ll«l) (B.37) Appendix B. Detailed Derivation of Pollaczek's Formula 154 Therefore, Ea = £ o e " | a | l / T. = E~0e-^ - juu.Je->a*fJL (e-l«-*/ll«l _ e-l-+-/ll«1) (B.38) (B.39) Using boundary condition (B.ll) EQ can be solved from the above two equations. Because dEa dy at = -\a\Eoe-M" = -\a\E0 y=0 t=0 y=0 (B.40) t=o (B.41) the following equation is got - —\a\E^ = - — (-\9\EZ - ju,ti.Ie-iam'e-t'W) Therefore, E0 is given by (B.42) ju>Ie-jax*e-ttW jule-jaxt e - t' \/a2+1/p2 EQ = (B.43) ±M + ±I'I i H + i>/5 + i7ri The final solution of 2?a and Ee can be got by applying the inverse Fourier transfor-mation to (B.38) and (B.39) with E0 being replaced by (B.43). Also, t and tf can now be changed back to — y and — y,, respectively. Ea and Ee are given by 1 r°° . Ea = —j E^da 27T J — O O /oo e — / v ei(x-xt)ada •°°iH + i \ / a 2 + 1/Pe2 / ======= cos((x — z,)aWa (B.44) Appendix B. Detailed Derivation of Pollaczek's Formula 155 1 r°° . Ee = — / Eee>*ada Jul r°° e - * / N / « 2 + i / r f - | « | t = -J— / — = ^x-x^ada 2TT 7-OO 2|«9| = / , = cos((x — x f)a)da / , e>lm-*Aadct (B.45) 2TT y_oo 2^/a2 + l/p 2 2?e in the above equation can be further simplified into juu,I ( r°° 2e(y+ytWa2+1/rt ^ •i-f- K0(D/pe) - MD'IVe) + I , cos((z - x,)a)da (B.46) Ee = where Ko(D/pe) = / e>lm—')ada (B.47) •>-«> 2yJo? + I/pi roo e-|-«-W/l\/ a 2 + 1/pJ , Ko(D7Pe) = / j Jl—'*da (B.48) D = y/(x-xfy + (y-y,)* (B.49) = yf(x-xf¥ + (y + yfy (B.50) K 0 is the zero order second kind modified Bessel funciton. Appendix C List of Symbols A, A : magnetic vector potential [A] : A vector Ab : A function values on FEM boundary [AB] : A vector of Dirichlet boundary nodes A B I : element of [AB] ak : distance to centre of PT cable of the A:th SC cable (&=A,B,C) an : unknown coefficient for trial function (pn O-mi • polynomial coefficient in Pm(Np,() K : value of A at global node n A(k) : A distribution in conductor k A(LCI) : A distribution in conductor k caused by conductor current J; [A^] : vector of A(ki) values in element Ei [AEI] = vector of A values at local nodes in element Ei element of [AEi] [Au] : A vector of unknown nodes B, B : magnetic field density [C] : shunt capacitance matrix per unit length of transmission lines C 0 i : coefficient in the formulas of SC coaxial cables C/,. : coefficient in the formulas of SC coaxial cables capacitance of the fcth insulation in a SC cable 156 Appendix C. List of Symbols GKi : coefficient in the formulas of SC coaxial cables Ct : diagonal element in [C] of SC coaxial cables C? : off diagonal element in [C] of SC coaxial cables C.i0 = direct self capacitance of conductor i Gmii : direct mutual capacitance between conductors i and j Cik : element of [C] Cp : perturbation coefficient of the earth d : differentiation [D*] : shape function differentiation matrix in element Ei D : charge density D : distance between a field point and a cable or diameter of strands D' : distance between a field point 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 Pm(Np, £) JtJ- : width of the jth. division from a surface of the ith conductor det[] : determinant of a matrix diag( ) : a diagonal matrix dl : integral element E, E : electrical field E0 : E value on the earth surface Ea : E field in the air Eb : E function values on FEM boundary Appendix C. List of Symbols EB : vector of node values of E\, in a FEM mesh Ec : E field in the earth with a deeply buried SC cable E. : E field in the earth EF : E field in the earth with a deeply buried current filament Ei : element i in a FEM mesh [Es] : source electrical field vector ESi : source electrical field in 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 Fmk ' elements in matrix [F] FB,„ : elements in matrix [FB] FEi • -1 mfe • Fmk in local node numbers in element Ei Fslh in local node numbers in element Ei [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 IeR, Iej : real and imaginary parts of the earth current I n : 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 5 hi '• real part and imaginary part of J; 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] : updates of [IE1] and [IE2] Iep : partial earth return current Iepc : partial earth return current of a deeply buried SC coaxial cable IePF : 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( tJ) values in element E\ Appendix C. List of Symbols 160 JR>JI '• * e a l and imaginary part of J Js : source current density [Js] : source current density vector Jsi : source current density in conductor i [JEi] : Jacobian transformation matrix in element Ei k : index K : number of conductors K0(), Ki() : second kind modified Bessel functions in 1st and 2nd order 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 NBB : number of Dirichlet boundary nodes in the earth : number of nodes in element Ei NP : 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 NP pe : earth complex penetration depth 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 rb : FEM boundary radius r e : inner earth radius or the outer radius of a stranded conductor req : 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 2,5 3 : areas of subtriangles Sc • cross-section area of a sector-shaped core conductor SR : solution region [Sc] ' conductor cross-section area matrix [S'c] • update of [Sc] [ScA] : a matrix modified from [Sc] Sck : cross-section area of conductor k [Scv] : a vector related to [EB] 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 Tmn : elements of [T] [TEi] : imaginary integral matrix for isoparametric element Ei : element of [TEi] [Ts] : imaginary integral matrix for simplex elements TSmn : elements of [Ts] [U] : coefficient matrix [UEi] : [U] in element Et Umn : elements of [U] Cfj, : Umn in local node numbers in element Ei Appendix C. List of Symbols 163 u z : unit vector of z axis in Cartician coordinate system 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 WEP ' electric energy found from the numerical field solution WEF : WEF with conductors i and j energized WEC : electric energy from circuit analysis with conductors i and j energized 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 xfiVi '• x a n d y of a buried filament [xEi] : vector of x coordinates of vertices for isoparametric element Ei xn,yn '• x and y coordinates of node n in a FEM mesh xp,yp : x and y coordinates of field point P yijjtejjfa : V coordinates of vertices in a simplex [yEi] : vector of y coordinates of vertices for isoparametric element 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 Ze : earth return impedance Zij : elements of [Z] Zki : submatrix in [Z] related to phases k and / (fc,/=A,B,C) Appendix C. List of Symbols 164 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)i : impedance related to the dielectrics between cylindrical conductors i and j [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 :mim2m3 '• shape function for simplex elements [/?] : shape function vector for isoparametric elements Pi : shape function for isoparametric elements S : real penetration depth in conductors 6i : S in conductor i 8e : 8 in the earth £() : Deric function 8 : real penetration depth in conductors e : permittivity €Ei '• permittivity in element Ei Cfc : permittivity of the k insulation cr : relative permittivity £ : 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 fie : permeability in the earth /*o : permeability in the vacuum fir : 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 : trial function for global node n VBi : trial function for Dirichlet node t : trial function for local node n in element Ei r : boundary surrounding the solution region r 0 : Dirichlet boundary Ti : homogeneous Neumann boundary rc> : 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 |
AggregatedSourceRepository | 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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0100550/manifest