J Control Autom Electr Syst (2017) 28:405–417 DOI 10.1007/s40313-017-0308-0 Modeling of Power Cables with Arbitrary Cross Section: From Parameter Calculation to Electromagnetic Transients Simulation Angelo Alfredo Hafner1 · Pablo Torrez Caballero2 · José Humberto A. Monteiro3 · Eduardo C. Marques Costa4 · Sérgio Kurokawa2 · Mauricio V. Ferreira Luz1 · Walter P. Carpes1 Received: 18 June 2016 / Revised: 19 September 2016 / Accepted: 12 February 2017 / Published online: 22 February 2017 © Brazilian Society for Automatics–SBA 2017 Abstract A full computer-based methodology is proposed for electromagnetic transient simulations in power cables characterized by an arbitrary cross-section geometry. The frequency-dependent parameters of the cables are calcu- lated using finite element method, and the three-phase cable modeling is carried out using modal decoupling and fitting techniques. The multiconductor representation of a sector- shaped cable is possible from the calculation of a constant and real modal transformation matrix, resulting four inde- pendent propagation modes (three phases and cable shield), which are modeled from the inclusion of frequency effect in the classic Bergeron method. The currents and voltages are expressed as a system of differential equations, which are presented as state equations and solved using numerical inte- B Angelo Alfredo Hafner angelo.hafner@posgrad.ufsc.br Pablo Torrez Caballero pablotorrezcaballero@gmail.com José Humberto A. Monteiro humbertomonteiro@ufac.br Eduardo C. Marques Costa educosta@pea.usp.br Sérgio Kurokawa kurokawa@dee.feis.unesp.br Mauricio V. Ferreira Luz mauricio.luz@ufsc.br Walter P. Carpes walter.carpes@ufsc.br 1 Universidade Federal de Santa Catarina – UFSC, Florianópolis, Brazil 2 Universidade Estadual Paulista – UNESP, São Paulo, Brazil 3 Universidade Federal do Acre – UFAC, Rio Branco, Brazil 4 Universidade de São Paulo – USP, São Paulo, Brazil gration methods. The proposed modeling technique allows the inclusion of time-variable and nonlinear elements during electromagnetic transient simulations in the time domain, which is not possible from frequency-domain models that are solved using inverse transforms. The proposed model is validated from results simulated using numerical Laplace transform and exact modal transformation matrix for calcu- lation of phase currents and voltages. Keywords Power cables · Electromagnetic transients · Frequency-dependent parameters · State equations 1 Introduction The impact of alternative power sources, distributed gen- eration and smart grids has been widely discussed in the present conjuncture of power systems in many countries. However, the proposal of alternative power sources (e.g., wind farms and solar energy) requires new technologies for integration of different power sources (distributed gen- eration) in the same conventional power grid by means of alternative techniques for power transmission and distribu- tion. For example, underground cables have been preferred to overhead transmission and distribution lines in several countries in Europe and North America due to environmen- tal and political issues (Kariyawasam et al. 2011; Orton and Samm 1997). Submarine and underground cables are widely used for interconnection of offshore wind farms and shore-based electrical distribution systems (Gouda and Dein 2015; Lazaropoulos 2014; Marshall et al. 2013; Kabalci et al. 2012; Papadopoulos et al. 2013). The development of other new offshore systems and renewable energy sources is being continuously encouraged as well as new techniques for electrical power transport using various configurations of 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s40313-017-0308-0&domain=pdf 406 J Control Autom Electr Syst (2017) 28:405–417 compact power cables with reduced cross section and non- conventional conductors’ shape. Low-voltage underground cables are widely used for power distribution in urban and populated areas in many countries. The power distribution by underground cables, instead overhead lines, provides several advantages in envi- ronmental and technical terms.Although underground power distribution iswell established inmany countries, a full atten- tion has been given for design of new power cables and impact of recent applications in smart grids on the distribu- tion grids using conventional underground cables (Gouda and Dein 2015; Lazaropoulos 2014;Marshall et al. 2013). An up- to-date research on new technologies, for underground power distribution, is the use of non-conventional compact cables for electrical power and telecommunication signal transmis- sion (Lazaropoulos 2014). It represents an asymmetric cable with four conductors: three sector-shaped conductors and a cylindrical shield conductor. The referred underground cable can provide a broadband communications platform for sev- eral applications in smart grids, such as advanced metering infrastructure, demand information, wide-area networking, storage/distribution grid management (Kabalci et al. 2012). Compact power cables are also applied for very specific situations, as for oil exploration on offshore platforms (Salles et al. 2010). Submarine machines are operated using umbili- cal cables composed of more than one independent three- phase power circuit with variable frequency and voltage. In addition, hydraulic steel tubes and independent conduc- tors for signal transmitting are integrated in the same power cable. Furthermore, submarine power cables are used for interconnection of offshorewind farms and shore-based elec- trical transmission/distribution systems (Chien and Bucknall 2006). Nowadays, sector-shaped cables have been used in low- and medium-voltage distribution grids because they present a smaller diameter for the same conductor cross section compared to cylindrical conductors of conventional coaxial and umbilical cables (Kariyawasam et al. 2011). However, sector-shaped cables could be used in several other appli- cations, such as interconnection of offshore wind farms and power/telecommunication transmission. An increasing effort is also observed in the development of new technologies applied for AC and DC power transmis- sion using underground coaxial cables and non-conventional sector-shaped cables. Modeling techniques and alternative methods for parameter calculation of power cables have been continuously included in the technical literature (Brito et al. 2016; Bonyadi-ram et al. 2012; Gudmundsdottir 2014). Non-conventional conductor shapes are designed in order to compact even more the cross section of power cables and also improving the power transmission capacity. On the other hand, power cables with reduced cross section and compact configurations require a more robust electrical insu- lation, alternative methodologies for parameter calculation and modeling techniques for electromagnetic transient sim- ulations (Bonyadi-ram et al. 2012). In this context, this paper proposes amulticonductor mod- eling using modal decoupling for electromagnetic transient simulations in sector-shaped cables. The model is devel- oped directly in the time domain using fitting techniques, without inverse transform and convolutions, which enables the inclusion of nonlinear and time-variable elements during simulations (switching, insulation failure, short circuit, etc.). The series impedances and admittance of the sector-shaped cable are calculated using the finite element method (FEM) because most analytical methods are restricted to conven- tional coaxial and umbilical cables (Kariyawasamet al. 2011; Bonyadi-ram et al. 2012; Hafner and Luz 2015). The pro- posedmethodology represents a full procedure for simulation of electromagnetic transient in power cables with arbitrary cross sections, since the parameter calculation to modeling using frequency-dependent equivalent circuits directly in the time domain. 2 Parameter Calculation of Sector-Shaped Cables using FEM Sector-shaped power cables are widely used for electri- cal power distribution in low- and medium-voltage levels. However, the non-conventional design of the cross section results in several restrictions in the analytical calculation of the impedance and admittance parameters. Thus, the finite element method (FEM) represents a powerful tool for cal- culation of power cables with non-conventional geometrical characteristics. The cross section of a sector-shaped cable is described in Fig. 1, and its geometrical and physical charac- teristics are presented in Table 1 (Ametani et al. 1980, 2015). Fig. 1 Sector-shaped cable 123 J Control Autom Electr Syst (2017) 28:405–417 407 Table 1 Geometrical and physical characteristics of the sector-shaped cable Component Radius (mm) σ (MS/m) ε (×ε0) μ (×μ0) Cores or phases r1 = 19.00 58.14 1.00 1.00 Inner insulator r2 = 25.00 0.000 4.10 1.00 Steel pipe r3 = 27.00 1.100 1.00 500 Outer insulator r4 = 30.00 0.000 2.30 1.00 Although the FEM represents an accurate procedure for parameter calculation of cables with arbitrary and non- conventional cross sections, there are well-known analytical methods in the technical literature. As a first example, there is an analytical method that proposes the approximation of an arbitrary cross-section conductor by an equivalent cylin- drical conductor using an analytical formulation to calculate the equivalent inner and outer radii (Ametani et al. 1980; Ametani and Fuse 1992). This analytical approach is not accurate for the entire frequency range considered in the pro- posed methodology, especially at high frequencies. The sub-conductor method represents another well-established analytical formulation for calculation of electrical parameters of non-conventional cables. This method was successfully applied for cylindrical conduc- tors and after adapted for sector-shaped cables. The cable is basically divided into three subsections (sector-shaped conductors) which are subdivided in the radial direction into elementary sub-conductors. This analytical method provides accurate results for frequencies no greater than 100 kHz, showing some inaccuracies at higher frequencies (Kariyawasam et al. 2011; Lucas and Talukdar 1978). On the other hand, a full-wave formulation for calculating the modal electromagnetic fields, and corresponding modal cur- rents and voltages, was recently proposed for multiconductor coaxial and sector-shaped cables. In contrast with the sub- conductor method, this last method shows a better accuracy for a wider range of frequencies (Habib and Kordi 2013). In fact, there are other important analytical methods in the technical literature for calculation of cable parameters with non-conventional cross sections, but only a few methods are accurate enough considering a wide range of frequencies. In this context, the FEM represents the most reliable procedure for validation of the proposed modeling technique. The electrical parameters of the sector-shaped cable in Fig. 1 are calculated using the FEM by means of two differ- ent computational programs: Gmsh and GetDP (Geuzaine and Remacle 2013; Dular and Geuzaine 2013). The Gmsh is a finite element mesh generator for pre- and post-processing applications with parametric input and advanced visualiza- tion tools (Geuzaine and Remacle 2013). The GetDP is an open-source finite element solver that can be applied for multi-physic problems (Dular andGeuzaine 2013). Thus, the Fig. 2 Cablemodeling (a) and sector-shaped conductormesh (b) using FEM first defines the mesh geometry, whereas the second defines the physical proprieties of the cable (insulation, conductors, etc.), boundary conditions and simulates the frequency- dependent current density andmagnetic flux in the phase and shield conductors. To calculate the series impedance, a mag- netic vector potential magnetodynamic formulation is used, and for the calculation of the parallel admittance, an elec- tric scalar potential electrostatic formulation is applied.More details about the procedure for calculation of the impedance and admittance parameters are described in “Appendix.” The first step to calculate the cable parameters using FEM is to model the sector-shaped cable in Fig. 1, based on geo- metric and physical characteristics provided in Table 1. The shape and mesh density are modeled using the Gmsh with a variable mesh density from the center of the sector-shaped conductor to the outer cross-section border, as described in Fig. 2. A frequency-dependent current density is observed through the cross sections of the sector-shaped conductors because of the skin and proximity effects (Brito et al. 2016; Ametani et al. 2015). Thus, a more detailed meshing is required closer to the cross-section borders of the sector- shaped conductors, as described in Fig. 2b. Since the sector-shaped cable is modeled using the Gmsh, the impedance and admittance parameters are calculated fol- lowing the numerical procedures described in “Appendix,” using the program GetDP. The self and mutual inductances are calculated from 10 Hz up to 1 MHz that represents an adequate frequency range for analysis of electromagnetic transients in power cables (Kariyawasam et al. 2011; Chien and Bucknall 2006; Ametani et al. 1980). Figure 3a shows a major current den- sity at the borders of the sector-shaped conductor because of the skin and proximity effects in power cables (Brito et al. 2016; Gudmundsdottir 2014). The global quantities like the voltage drop or the current intensity, coupling the cable with an external electrical cir- cuit, are often those which permit to describe the overall 123 408 J Control Autom Electr Syst (2017) 28:405–417 Fig. 3 Cablemodeling (a) and sector-shaped conductormesh (b) using FEM Fig. 4 Self and mutual resistances of the sector-shaped cable functioning of the system. In this paper, the frequency- dependent resistances and inductances are calculated using a circuit relation associated with each conductor (cores and armor). Figures 4 and 5 show the resistances and the induc- tances obtained by FEM, respectively. In Fig. 4, the curve raa is the self-resistance of the phase a (impressed current), rab is the mutual resistance between a and b, whereas rag is the mutual resistance between a and Fig. 5 Self and mutual inductances of the sector-shaped cable Fig. 6 Electric field mapping obtained using the FEM armor g. The self-resistance of the armor rgg is also described in Fig. 6. It is important to emphasize that rab is similar to rac and rbc because the sector-shaped cable has a symmetrical cross section. The same statement is valid for self and mutual inductances in Fig. 5. The capacitance matrix is calculated based on the relation between electric charge (C/m) and potential difference (V ) impressed between the elements, as properly described in “Appendix.” Figure 6 shows that the electric field is not homogeneously distributed in the cross section of the conductors. The electric field intensity is more pronounced at the three corners of the cross section. Therefore, from the electric field calculated based on the circuit diagram in “Appendix,” the capacitance matrix can be calculated following the well-established analytical for- mulation in the technical literature (Hafner and Luz 2015; Ametani et al. 2015): 123 J Control Autom Electr Syst (2017) 28:405–417 409 [C] = ⎡ ⎢⎢⎣ 601.35 −164.07 −164.07 −273.20 −164.07 601.35 −164.07 −273.20 −164.07 −164.07 601.35 −273.20 −273.20 −273.20 −273.20 2, 034.18 ⎤ ⎥⎥⎦ (1) The self and mutual capacitances in (1) present minor vari- ations compared to capacitance matrix obtained using the sub-conductor method for the same sector-shaped cable in reference (Kariyawasam et al. 2011). The capacitance mag- nitudes are expressed in ηF/km. 3 Multiconductor Representation Using Modal Decoupling The modal decoupling represents a well-established tech- nique for modeling of multiconductor overhead transmission lines and underground cables. This technique has been widely used for representation of three-phase power systems, which are decoupled into three independent propagation modes. This means that each propagation mode can be mod- eled as a single-phase line in the frequency domain as a two-port circuit or directly in the time domain using fitting techniques (Costa et al. 2010, 2013; Caballero et al. 2015). The three-phase system decoupling is the most usual method to model symmetrical and untransposed transmis- sion systems, which can be solved without the explicit and complex representation of mutual parameters and electro- magnetic coupling between conductors (Costa et al. 2013). The modal decoupling of the impedance and admittance matrices of the system is carried out by a modal transfor- mation matrix, as described in (2) and (3) (Costa et al. 2013). [Zm] = [T ]T [Z ] [T ] (2) [Ym] = [T ]T [Z ] [T ] (3) Terms [Z] and [Y ] are impedance and admittance matrices of a multiconductor line or cable. The modal matrix is denoted as [T ], where [T ]T and [T ]−1 are the transposed and inverse forms, respectively. The columns of the modal transforma- tion matrix [T ] are eigenvectors obtained from eigenvalues of the matrix product [Y ][Z], which are calculated using the well-known Newton–Raphson method (Costa et al. 2013). Applying (2) and (3) to decouple the sector-shaped cable in Sect. 2, four propagationmodes are obtained and thematri- ces [Zm] and [Ym] can be expressed as: [Zm] = ⎡ ⎢⎢⎣ Z1 0 0 0 0 Z2 0 0 0 0 Z3 0 0 0 0 Z4 ⎤ ⎥⎥⎦ ; [Ym] = ⎡ ⎢⎢⎣ Y1 0 0 0 0 Y2 0 0 0 0 Y3 0 0 0 0 Y4 ⎤ ⎥⎥⎦ (4) From (4), it is possible to verify that each propagation mode is completely decoupled from others because there are no mutual terms in the modal impedance and admittance matri- ces [Zm] and [Ym]. Thus, each propagation model can be represented as a single-phase line using various frequency- or time-domain models in the technical literature for simulation of electromagnetic transients (Costa et al. 2010; Caballero et al. 2015; Gómes and Uribe 2009). If each propagation mode is approached as a single-phase line, the modal currents and voltages (modal domain) can be expressed in vector form as follows: [Im] = [ I1 I2 I3 I4 ] (5) [Vm] = [ V1 V2 V3 V4 ] (6) There are basically two procedures to calculate the currents and voltages in (5) and (6), respectively. The first is based on the two-port representation of each propagationmode, which the currents and voltages are calculated in the frequency domain and converted to the time domain using inverse trans- forms and convolutions (Gómes andUribe 2009). The second procedure is carried out directly in the time domain, using fitting techniques and numerical integration methods (Costa et al. 2010; Caballero et al. 2015). The modeling directly in the time domain presents several advantages compared to the method by two-port representation, such as the inclusion of nonlinear and time-variable elements during simulations. In sequence, after calculation of the currents and voltages at the four propagation modes of the sector-shaped cable, these values are converted to the phase domain using amodal transformationmatrix, as expressed in (7) (Costa et al. 2013). [ Iph ] = [T ] [Im] T ; [ Vph ] = [T ]−T [Vm] T (7) The vectors [Iph] and [Vph] are the currents and voltages at the cores a, b, c and at the armor g of the sector-shaped cable, as expressed in (8) and (9). [ Iph ] = [ Ia Ib Ic Ig ] (8)[ Vph ] = [ Va Vb Vc Vg ] (9) The modal transformation in (7) can be accomplished in frequency-domain models using the same transformation matrix [T ], which is composed of complex and frequency- dependent values. However, the same modal transformation matrix cannot be used for line/cable models in the time domain. In this case, the real and constant Clarke’s matrix is an alternative approach for modeling of three-phase systems, which is not the case of the proposed sector-shaped cable that presents four propagation modes (Costa et al. 2010, 2013). 123 410 J Control Autom Electr Syst (2017) 28:405–417 Fig. 7 Impedance values of each propagation mode of the sector- shaped cable: modal resistance Rm (a) and modal reactance Xm (b) 4 Fitting a Modal Transformation Matrix for Cable Modeling in the Time Domain The modal matrices in (4), [Zm] and [Ym] can be obtained as a function of [Z] and [Y ], using the exact transformation matrix [T ] that was introduced in Sect. 3 in (2) and (3). The modal resistance and reactance that compose the impedance of each propagationmode of the sector-shaped cable in Fig. 1 are described in Fig. 7. Figure 7 shows the resistance and reactance values of the modal impedances that compose the matrix [Zm], with four propagation modes resulted from the modal decoupling of the sector-shaped cable in Fig. 1. The real and imaginary terms of the modes 1 and 2 are practically overlapped, as shown in Fig. 7a, b, respectively. The modal admittance matrix [Ym] is composed only by imaginary terms that represent the capacitive reactances of each propagation mode, where the modal capacitances are constant, such as expressed in the phase domain in (1). This way, [Ym] can be expressed by a constant capacitance matrix in ηF/km: [Cm] = ⎡ ⎢⎢⎣ 765.41 0 0 0 0 765.41 0 0 0 0 1210.10 0 0 0 0 1097.40 ⎤ ⎥⎥⎦ (10) The modal admittance matrix can be expressed as a function of the capacitance matrix [Cm] and the angular velocity ω as follows: Fig. 8 Real components of the exact modal transformation matrix [T ] [Ym] = jω [Cm] (11) All propagation characteristics (propagation function, wave propagation velocity and impedance characteristic of each propagation mode) can be obtained from the impedance and admittance matrices [Zm] and [Ym]. The four propagation modes of the cable can be obtained using the matrix [T ], but the conversion of the modal cur- rents and voltages to phase values in (7) cannot be carried out from a transformation matrix composed of complex and frequency-dependent terms, because the cable model- ing is developed directly in the time domain and requires a real and constant transformation matrix with dimension four by four. In this context, the behavior of the complex and frequency-dependent terms in [T ], obtained from the impedance and admittance parameters of the sector-shaped cable, is evaluated in order to find a suitable real and constant transformation matrix [Tc]. The real parts of each compo- nent of [T ], obtained using the Newton–Raphson method, are described as a function of the frequency in Fig. 8. The frequency range in Fig. 8 is from10Hz up to 0.3MHz, for a more detailed observation of the frequency-dependent behavior of the real components in [T ]. These components are practically constant for frequencies up to10Hzandhigher than 0.3MHz,which are not explicit in the frequencywindow in Figs. 7 and 8 in order to highlight the variations between 200 Hz and 2kHz. Figure 8 shows that some real components present oscilla- tions in short-frequency intervals between 100Hz and 1 kHz. In sequence, for frequencies higher than 1 kHz, all real com- 123 J Control Autom Electr Syst (2017) 28:405–417 411 Fig. 9 Imaginary components of the exact modal transformation matrix [T ] ponents converge to constant values in the great frequency range up to 1 MHz. On the other hand, Fig. 9 shows that the imaginary components of the first two rows oscillate around zero from 10 Hz up to 0.3 MHz. For frequencies higher than 10 kHz, all imaginary components in [T ] are approximately null up to 1 MHz. The imaginary components of the last two rows in [T ] showed the same low oscillations and converging to zero, as in Fig. 9. From these observations, constant and real values can be chosen to compose an approximated transformation matrix, whereas the imaginary components can be neglected. Another possibility is to fit a constant and real modal trans- formationmatrix based on the electromagnetic transient to be simulated. For example, for transient simulations composed of low frequencies (e.g., switching impulse), real compo- nents between 10 and 200 Hz are selected from Fig. 8 to compose an approximated transformation matrix [Tc] with real and constant terms, as proposed in (12). [Tc] = ⎡ ⎢⎢⎣ 0.79 0.19 0.59 0.05 −0.57 0.59 0.59 0.05 −0.22 0.78 0.59 0.05 0 0 0.23 1.05 ⎤ ⎥⎥⎦ (12) For transient simulations characterized by a wide range of frequencies, such as a steep-front impulse, real component values above 1 kHz should be considered (Fig. 7), which can be also assumed as constant terms in the frequency range from 1 kHz up to 1 MHz. Thus, for simulation of fast and impulsive transients (e.g., atmospheric impulse), the approx- imated modal transformation matrix [Tc] can be fitted using high-frequency terms from figure, as expressed in (13). [Tc] = ⎡ ⎢⎢⎣ 0.75 0.49 1.0 1.15 −0.55 0.47 1.0 1.15 −1.15 0.9 1.0 0.05 0 1.0 1.0 0 ⎤ ⎥⎥⎦ (13) From the proposed approach for the modal transformation matrix, (7) can be expressed as: [ Iph ] = [Tc] [Im] T ; [ Vph ] = [Tc] −T [Vm] T (14) The proposed cable model and methodology for calcula- tion of an approximated transformation matrix [Tc] can be validated by comparison with results obtained using the exact frequency-dependent matrix [T ], cablemodeling based on the two-port representation of each propagation mode and numerical Laplace transform (NLT) (Gómes and Uribe 2009). 5 Frequency-Dependent Modeling in the Time Domain The frequency-dependent modeling of the sector-shaped cable is carried out using modal decoupling in the multi- conductor representation, as described in Sect. 4, and fitting techniques for inclusion of the frequency effect in the well- known Bergeron method. The original Bergeron method represents the active losses of the system by lumped resistances located at the send- ing and receiving ends of a single electrical circuit with a historical current source. In the proposed Bergeron circuit, these lumped resistances are substituted by an equivalent cir- cuit that is modeled from a rational function that fits the frequency-dependent parameters of each propagation mode of the cable. The distributed parameters of each mode can be represented as a cascade of frequency-dependent Berg- eron circuits and state equations, which are solved using any numerical integration method. From this procedure, the sector-shaped cable can be modeled directly in the time domain taken into account nonlinear and time-variable ele- ments during electromagnetic transient simulations. The proposed time-domain model can be validated based on sim- ulations obtained from the cable modeling by equivalent two-port circuits and NLT. 5.1 Frequency-Dependent Bergeron Circuit The cable modeling is carried out from the frequency- dependent parameters fitting in the equivalent Bergeron circuit, as described in Fig. 10. Thus, each propagation mode of the cable is modeled as a single-phase line as a function of the frequency (Caballero et al. 2015). Terms Ik and Im are equivalent current sources at the sending and receiving ends of the equivalent circuit, respec- tively. Sources Ik and Im are known at the state time t from the past history at time (t − τ ). Terms vk and vm are the voltages at the sending and receiving ends of the equivalent circuits, respectively. The impedance characteristic of each propagation mode Z0 is expressed as a function of the shunt capacitance and inductance in direct current. The RL circuits atm and n are obtained froma rational function that is applied 123 412 J Control Autom Electr Syst (2017) 28:405–417 Fig. 10 Frequency-dependent Bergeron circuit for fitting the series impedance of the cable, as expressed in (15) (Gustavsen and Semlyen 1998). Z (ω) ≈ Zfit (ω) = R′ 0 + jωL ′ 0 + ∑n i=1 ( jωR′ i( jω + R′ i/L ′ i ) ) (15) Termω is the angular velocity. The per-unit-length resistance R′ 0 and per-unit-length inductance L ′ 0 are values for ω = 0. The frequency range is fitted as a function of the quantity of parallel RL circuits in the frequency-dependent Bergeron circuit in Fig. 10. There are several fitting procedures available in the tech- nical literature to obtain an approximated rational function from a general function, such as Z(ω). However, the vector fitting algorithm shows to be accurate and robust for smooth and resonant responses with high-order and wide frequency bands (Gustavsen and Semlyen 1998). Considering the frequency-dependent Bergeron circuit in Fig. 10, indicated by the impedance Z(ω)/2, the relationship of voltage on the resistors and inductors at the node k is expressed in (16)–(18): R′ 1 2 (ik0 − ik1) = L ′ 1 2 dik1 dt (16) R′ 2 2 (ik0 − ik2) = L ′ 2 2 dik2 dt (17) R′ n 2 (ik0 − ikn) = L ′ n 2 dikn dt (18) In (16)–(18), terms ik0, ik1 to ikn are the currents in the resis- tor R′ 0/2 and inductors L ′ 1/2 to L ′ n /2, respectively. Following the Kirchoff’s current law, the current through the resistors R′ 1/2 and R′ n /2 are (ik0 – ik1) and (ik0 – ikn), respectively. Thus, the general expression is developed associating the past-history current source Ik(t − τ) with the state current ik0, which is equivalent to the current ik,m(t) in Fig. 10. Vk − R0ik0 − R1 (ik0 − ik1) − · · · − Rn (ik0 − ikn) +Z0 Ik (t − τ) = Z0ik0 (19) The expressions of the two branches are similar, considering that the impedance Zfit is divided by two, as described in Fig. 10. In order to simplify the formulation in (19), the per-unit- length parameters R′ 0/2 and L ′ 0/2 are denoted as R0 and L0. The same notation is considered for n parallel RL circuits. 5.2 Representation by State Equations Figure 10 shows that the frequency-dependent Bergeron circuit is composed of two similar circuits represented by (16)–(19), which can be both represented as state equations: [ . Ik ] = [Ak] [Ik] + [Bk] [Sk] (20) Vector [Ik] represents the currents in L1 to Ln : [Ik] = [ik1 ik2 . . . ikn] T (21) Thus, the derivative of the currents in [Ik] can be expressed: [ . Ik ] = [ ik1 dt ik2 dt . . . ikn dt ] (22) The state matrix [Ak] is expressed in (23) that is a function of the characteristic impedance Z0 of the Bergeron circuit and the lumped elements in Fig. 10. [Ak] = ⎡ ⎢⎢⎢⎢⎢⎢⎣ R1 L1 ( −1 + R1 Z0+∑n i=1 Ri ) R1 L1 ( R2 Z0+∑n i=1 Ri ) · · · R1 L1 ( Rn Z0+∑n i=1 Ri ) R2 L2 ( R1 Z0+∑n i=1 Ri ) R2 L2 ( −1 + R2 Z0+∑n i=1 Ri ) · · · R2 L2 ( Rn Z0+∑n i=1 Ri ) ... ... . . . ... Rn Ln ( R1 Z0+∑n i=1 Ri ) Rn Ln ( R2 Z0+∑n i=1 Ri ) · · · Rn Ln ( −1 + Rn Z0+∑n i=1 Ri ) ⎤ ⎥⎥⎥⎥⎥⎥⎦ (23) 123 J Control Autom Electr Syst (2017) 28:405–417 413 Analogously as in (23), the matrix [Bk] is expressed based on Z0 and the RL circuit in Fig. 10: [Bk] = ⎡ ⎢⎢⎢⎢⎢⎢⎣ R1 L1 ( 1 Z0+∑n i=1 Ri ) R1 L1 ( Z0 Z0+∑n i=1 Ri ) R2 L2 ( 1 Z0+∑n i=1 Ri ) R2 L2 ( Z0 Z0+∑n i=1 Ri ) ... ... Rn Ln ( 1 Z0+∑n i=1 Ri ) Rn Ln ( Z0 Z0+∑n i=1 Ri ) ⎤ ⎥⎥⎥⎥⎥⎥⎦ (24) Vector [Sk] is composed of the voltage source Vk and the historical current Ik(t– τ): [Sk] = [ Vin (t) Ik (t − τ) ] (25) Term Vin(t) is the input voltage source applied at the sending end of the Bergeron circuit. From the state-space formula- tion in (22), ik,m(t) and Vk(t) can be analytically expressed from elements in [Ik] and [Sk], as expressed in (26) and (27), respectively (Caballero et al. 2015). ik,m (t) = ∑n q=1 ( Rqikq Z0 + ∑n i=1 Ri ) + 1 Z0 + ∑n i=1 Ri Vk + Z Z0 + ∑n i=1 Ri Ik (t − τ) (26) Vk (t) = Vin (t) − ik,m (t) ∑n i=1 Ri (1 + iki ) (27) The formulation of the branch m is similar to k, as described in Fig. 10. However, if an impedance ZL is connected to the receiving end at m, the state matrix [Am] can be expressed as: [Am] = ⎡ ⎢⎢⎢⎢⎢⎢⎣ R1 L1 ( −1 + R1 Z0+ZL+∑n i=1 Ri ) R1 L1 ( R2 Z0+ZL+∑n i=1 Ri ) · · · R1 L1 ( Rn Z0+ZL+∑n i=1 Ri ) R2 L2 ( R1 Z0+ZL+∑n i=1 Ri ) R2 L2 ( −1 + R2 Z0+ZL+∑n i=1 Ri ) · · · R2 L2 ( Rn Z0+ZL+∑n i=1 Ri ) ... ... . . . ... Rn Ln ( R1 Z0+ZL+∑n i=1 Ri ) Rn Ln ( R2 Z0+ZL+∑n i=1 Ri ) · · · Rn Ln ( −1 + Rn Z0+ZL+∑n i=1 Ri ) ⎤ ⎥⎥⎥⎥⎥⎥⎦ (28) The currents through L1 to Ln of both k andm are similar, which means that [Im] and [Ik] are: [Im] = [Ik] = [ik1 ik2 . . . ikn] T (29) The matrix [Bm] is a single-column vector, since the node m is not connected to a voltage source, where the impedance ZL represents an equivalent Thevenin circuit: [Bm] = ⎡ ⎢⎢⎢⎢⎢⎢⎣ R1 L1 ( Z0 Z0+ZL+∑n i=1 Ri ) R2 L2 ( Z0 Z0+ZL+∑n i=1 Ri ) ... Rm Lm ( Z0 Z0+ZL+∑n i=1 Ri ) ⎤ ⎥⎥⎥⎥⎥⎥⎦ (30) Thematrix [Sm] contains the history current Im(t– τ), voltage Vm and current im,k , as in (31). Vm (t) = Vout (t) − im,k (t) ∑n i=1 Ri (1 + iki ) (31) ik,m (t) = ∑n q=1 ( Rqikq Z0 + ZL + ∑n i=1 Ri ) + Z0 Z0 + ZL + ∑n i=1 Ri Im (t − τ) (32) The voltage on the impedance ZL is represented by the output source Vout(t), expressed in (33). Vout (t) = −im,k (t) ZL (33) The state matrices in (20), for both nodes k and m in Fig. 10, are dependent of the fitted R and L parameters and the char- acteristic impedance Z0. The state equations can be solved using integration methods (Costa et al. 2010; Caballero et al. 2015). 5.3 Cascade of Frequency-Dependent Bergeron Circuits The line or cable approach by cascade of equivalent circuits can represent accurately the distributed characteristic of the electrical parameters. Furthermore, the frequency response of the model is proportional to the number of equivalent cir- cuits used in cascade (Costa et al. 2010). Thus, the cascade representation shows to be a good approach formodeling and simulation of fast and impulsive electromagnetic transients, which are composed of a wide range of frequencies. In Sect. 5.2, the state equations represent a single frequency-dependent Bergeron circuit. In this section, the state matrices are introduced for a cascade with h equivalent frequency-dependent Bergeron circuits, as in Fig. 11. Each block in Fig. 11 represents a frequency-dependent Bergeron circuit (Fig. 10). The cascade is represented by (h+ 123 414 J Control Autom Electr Syst (2017) 28:405–417 Fig. 11 Cascade of frequency-dependent Bergeron circuits 1) branches, where the input branch connected to Vin is Zfit/2 and the output branch in the equivalent circuit h is also Zfit/2. The intermediary interconnected branches from circuit 1 to h represent Zfit. For example, the receiving end m of the equivalent Bergeron circuit 1 is Zfit/2 and the sending end k of the circuit 2 is also Zfit/2, and the interconnection of the equivalent circuits 1 and 2 results in Zfit. Thus, the current im,k from the receiving end of (h– 1) is similar to ik,m at the sending end of h: i (h) k,m = −i (h−1) m,k (34) The connection of the RL circuits in (h– 1) and h results in some variations in the state matrix [B] and vector [S], which are renamed as [Bmk] and [Smk] in (35) and (36), respectively. The state matrix [A] remains as [Ak]. [Bmk] = ⎡ ⎢⎢⎢⎢⎢⎢⎣ R1 L1 ( Z0 2(Z0+∑n i=1 Ri) ) R1 L1 ( −Z0 2(Z0+∑n i=1 Ri) ) R2 L2 ( Z0 2(Z0+∑n i=1 Ri) ) R2 L2 ( −Z0 2(Z0+∑n i=1 Ri) ) ... ... Rn Ln ( Z0 2(Z0+∑n i=1 Ri) ) Rn Ln ( −Z0 2(Z0+∑n i=1 Ri) ) ⎤ ⎥⎥⎥⎥⎥⎥⎦ (35) [Smk] = [ I (h−1) m (t − τ) I (h) k (t − τ) ] (36) Thevector [Smk] is composedof the historical current sources Im(t– τ) in the Bergeron circuit (h– 1) and Ik(t– τ) in the circuit h, as subscripted in the terms in (36). From (35) and (36), (20) can be reformulated for two con- nected Bergeron circuits: [ İk ]h−1 = [Ak] [Ik] h−1 + [Bmk] [Smk] (37) Analogously to (26) and (31), the current flowing from h– 1 to h can be analytically expressed: i (h−1) m,k = −i (h) k,m = ∑n q=1 ( Rqikq Z0 + ∑n i=1 Ri ) + Z0 2 ( Z0 + ∑n i=1 Ri ) ( I (h−1) m (t − τ) − I (h) m (t − τ) ) (38) The voltages Vm(t) at the circuit h– 1 and Vk(t) at h are expressed as: V (h−1) m (t) = Z [ ik0 − I (h−1) m (t − τ) ] (39) V (h) k (t) = −Z [ ik0 + I (h) k (t − τ) ] (40) For the input branch, [Ak], [Bk] and [S], in (23), (24) and (25), respectively, are considered as in (20). For the output branch, the state matrices [Am] and [Bm] are expressed in (28) and (30), respectively. At the last Bergeron circuit of the cascade, [S] is represented solely by the history current Im(t– τ). Based on themodeling technique described in this section, eachpropagationmodeof the sector-shaped cable (Fig. 1) can be modeled as a cascade of frequency-dependent Bergeron circuits. The transient currents and voltages for each prop- agation modes are calculated and then converted to phase values using (14) and the constant and real transformation matrices in (12) and (13). 6 Time-Domain Simulations Using the Proposed Model The electromagnetic transient simulations obtained from the proposed modeling methodology are compared to simula- tions using NLT and the exact frequency-dependent trans- formation matrix [T ]. The first principal difference between these two modeling methodologies is that the propagation modes of the proposed modeling technique are represented directly in the time domain by fitting techniques, using a cascade of frequency-dependent Bergeron circuits, whereas the methodology using NLT represents each mode as a two-port circuit in the frequency domain. The second dif- ference between the modeling methodologies is that the phase currents and voltages simulated from the proposed methodology are obtained from the modal domain using an approximated real and constant transformation matrix [Tc], (12)–(14), whereas the same values in the modeling by NLT are obtained using the exact modal matrix [T ] in (7). The proposed modeling methodology can be validated compar- ing simulations obtained using the NLT, as reference results, and simulations carried out by the proposed model. Figure 12 describes the sending end of the phase/core a of the sector-shaped cable connected to a voltage source U (t) and other two phases and armor g connected at the sending end to the ground. The three phases and armor are open at the receiving end. The sector-shaped cable, represented in Fig. 12, is 100 km long and has the same characteristics and parameters calculated in Sect. 2. The first simulation consists of an unitary voltage step applied at the sending end of the phase a, by means of the 123 J Control Autom Electr Syst (2017) 28:405–417 415 Fig. 12 Validation of the proposed modeling methodology Fig. 13 Voltage transients at the receiving end of the phases a, b, c and armor g simulated from a switching impulse: NLT (dotted curves) and proposed modeling (solid curves) voltage source U (t). Results obtained from the proposed model, which each propagation mode is represented by 100 Bergeron circuits in cascade, are simulated using the modal transformation matrix in (12) that is valid for frequencies up to approximately 400–500 Hz. The voltages at the receiv- ing end of the sector-shaped cable, simulated using the two models, are described in Fig. 13. The dotted curves represent the results obtained from the line model using NLT, and the solid curves are the transient voltages simulated using the proposed time-domain model- ing. The phases a, b, c and the cable armor g are explicitly indicated in Fig. 13. From the same configuration in Fig. 12, the transient voltages at the receiving end of the cable are simulated con- sidering an atmospheric impulse applied at the sending end of the phase a (1.2/50 μs voltage impulse). The modal trans- formation matrix in (13) is applied for calculation of the transient voltages at the receiving end of the cable. Dotted and solid curves are results obtained using NLT and the pro- posed model, respectively, as described in Fig. 14. Figures 13 and 14 show that simulations obtained from the proposed model, using fitting techniques for calculation of a real and constant modal transformation matrix and a frequency-dependent equivalent circuit of each propagation mode, are almost similar to the reference values obtained using the cable model based on the exact transformation matrix and NLT. Fig. 14 Voltage transients at the receiving end of the phases a, b, c and armor g simulated from an atmospheric impulse: NLT (dotted curves) and proposed modeling (solid curves) 7 Conclusions A modeling methodology for power cables with arbitrary cross sections was proposed for electromagnetic transient simulations. The proposed methodology presents a detailed procedure and discussions that approach since the parameter calculation of sector-shaped cables to the time-domain tran- sient simulations using well-established and novel modeling techniques, such as the calculation of a real and constant modal transformation matrix for cables with more than three propagation modes and representation using cascade of frequency-dependent Bergeron circuits. The proposed model was validated based on results simu- lated using NLT and the exact modal transformation matrix for calculation of the transient currents and voltages. The cable representation using an approximated modal transfor- mation matrix and frequency-dependent Bergeron circuits proved to be accurate and without numerical oscillations. In addition, the proposed modeling methodology can be extended for conventional umbilical and coaxial power cables with more than three propagation modes, without restrictions in the shape of the transversal cross section of the phase conductors. Acknowledgements São Paulo Research Foundation – FAPESP (Procs. 14/17051-0 and 15/10204-8) and National Counsel of Tech- nological and Scientific Development – CNPq (306142/2015-5) Appendix: Impedance and Admittance Parameters The numerical approach for calculation of the self andmutual parameters of multiconductor cables, including phase con- ductors and shieldmesh, can be carried out applying a current of 1.0 A at one of the conductors and calculating the volt- age at all conductors. In the case of a sector-shaped cable, there are four conductors: three-phase cores and shield wire (armor). The self-impedances are calculated dividing the induced voltage by the impressed current,whereas themutual 123 416 J Control Autom Electr Syst (2017) 28:405–417 Fig. 15 Cable representation for calculation of the impedance param- eters are calculated dividing the induced voltage on other elements without impressed current (Hafner and Luz 2015). This pro- cedure is illustrated in Fig. 15. Initially, considering a current of 1.0 A through the phase conductor a, the mutual parameters can be calculated between a and other two phase conductors and armor. The procedure is valid considering the same current imposed in phases b or c. Thus, the frequency-dependent impedances of the sector-shaped cable can be calculated using the FEM, as a function of the current density (A/m2) and the magnetic flux (Wb/m). Figure 15 shows that the cable is composed of phase cores a, b and c and armor wire g. Thus, from the voltage drops along the sector-shaped cable, the impedance matrix can be expressed in (41) (Ametani et al. 2015). ⎡ ⎢⎢⎣ Vag − Va′g′ Vbg − Vb′g′ Vcg − Vc′g′ Vgg − Vg′g′ ⎤ ⎥⎥⎦ = ⎡ ⎢⎢⎣ Zaa Zab Zac Zag Zba Zbb Zbc Zbg Zca Zcb Zcc Zcg Zga Zgb Zgc Zgg ⎤ ⎥⎥⎦ ⎡ ⎢⎢⎣ Ia Ib Ic Ig ⎤ ⎥⎥⎦ (41) Terms Zaa, Zbb, Zcc are the self-impedances of the phase conductors,whereas Zab, Zac and Zbc aremutual impedances. The impedances Zag, Zbg and Zcg are the mutual terms between phases a, b and c and the shield/ground wire g, respectively. The self and mutual capacitances are considered as invari- able in power cables and overhead transmissions lines up to very high frequencies. Therefore, the capacitances are usually considered constant in transmission line and cable modeling for electromagnetic transient analysis (Ametani et al. 2015). The self and mutual capacitances are calculated as a func- tion of the electric charge and the voltage at the conductors. Thus, the capacitance matrix of the sector-shaped cable can be expressed as (Hafner and Luz 2015; Ametani et al. 2015): Fig. 16 Circuit representation for calculation of the capacitances ⎡ ⎢⎢⎣ Ia − Ia′ Ib − Ib′ Ic − Ic′ Ig − Ig′ ⎤ ⎥⎥⎦ = ⎡ ⎢⎢⎣ Caa Cab Cac Cag Cba Cbb Cbc Cbg Cca Ccb Ccc Ccg Cga Cgb Cgc Cgg ⎤ ⎥⎥⎦ ⎡ ⎢⎢⎣ Vag Vbg Vcg Vgg ⎤ ⎥⎥⎦ (42) The voltage and current vectors at the left hand in (41) and (42) are implicitly represented in terms of current density, magnetic flux and electric field, which are obtained using FEM and following the proposed boundary conditions for impedance and capacitance calculation in Figs. 15 and 16, respectively. References Ametani, A. (1980). A general formulation of impedance and admit- tance of cables, IEEE Trans. Power Apparatus and Systems, 3, PAS-99, pp. 902–910. Ametani, A., & Fuse, I. (1992). Approximate method for calculating the impedances of multiconductors with cross sections of arbitrary shapes. Electrical Engineering in Japan, 112(2), 117–123. Ametani, A., Ohno, T., & Nagaoka, N. (2015).Cable system transients: Theory, modeling and simulation. New York: Wiley-IEEE Press. Bonyadi-ram, S., Kordi, B., & Bridges, G. E. (2012). A full-space conformalmapping for the calculation of series impedance of over- head transmission lines and underground cables. Electric Power System Research, 91, 95–103. Brito, A. I., Machado, V. M., Almeida, M. E., & Neves, M. G. (2016). Skin and proximity effects in the series-impedance of three-phase underground cables. Electric Power Systems Research, 130, 132– 138. Caballero, P. T., Costa, E. C. M., & Kurokawa, S. (2015). Frequency- dependent multiconductor line model based on the Bergeron method. Electric Power System Research, 127, 314–322. Chien, C. H., & Bucknall, R. W. G. (2006). Theoretical aspects of the harmonic performance of subsea AC transmission systems for offshore power generation schemes. IEE Proceedings Generation, Transmission and Distribution, 153(5), 599–609. Costa, E. C. M., Kurokawa, S., Pinto, A. J. G., Kordi, B., & Pisso- lato, J. (2013). Simplified computational routine to correct modal decoupling in transmission lines and power systemsmodeling. IET Science Measurement and Technology, 7(1), 7–15. 123 J Control Autom Electr Syst (2017) 28:405–417 417 Costa, E. C. M., Kurokawa, S., Pissolato, J., & Prado, A. J. (2010). Efficient procedure to evaluate electromagnetic transients on three-phase transmission lines. IET Generation Transmission and Distribution, 4, 1069–1081. Dular, P., & Geuzaine, C. (2013). GetDP Reference Manual. Geuzaine, C., & Remacle, J. (2013). Gmsh Reference Manual. Gómes, P., & Uribe, F. A. (2009). The numerical Laplace transform: An accurate technique for analyzing electromagnetic transients on power system devices. International Journal of Electrical Power and Energy Systems, 31(2–3), 116–123. Gouda, O. E., & Dein, A. Z. E. (2015). Improving underground power distribution capacity using artificial backfill materials. IET Gener- ation Transmission and Distribution, 9(15), 2180–2187. Gudmundsdottir, U. S. (2014). Proximity effect in fast transient simula- tions of an underground transmission cable.Electric Power System Research, 115, 50–56. Gustavsen, B., & Semlyen, A. (1998). Combined phase and modal cal- culation of transmission line transients based on vector fitting. IEEE Transactions on Power Delivery, 13(2), 596–604. Habib, S., & Kordi, B. (2013). Calculation of multiconductor under- ground cables high-frequency per-unit-length parameters using electromagnetic modal analysis. IEEE Transactions on Power Delivery, 28(1), 276–284. Hafner, A. A., &Luz,M.V. F. (2015). Impedance and admittance calcu- lations of a three-core power cable by the finite element method. In International conference on power systems transients—IPST, Cavtat, Croatia. Kabalci, E., Kabalci, Y., & Develi, I. (2012). Modelling and analysis of a power line communication system with QPSK modem for renewable smart grids. International Journal of Electrical Power and Energy Systems, 34(1), 19–28. Kariyawasam, K. K. M. A., Gole, A. M., & Kordi, B. (2011). Accu- rate electromagnetic transient modeling of sector-shaped cables. In International conference on power systems transients—IPST, Delft, Netherlands. Lazaropoulos, A. G. (2014). Numerical evaluation of broadband trans- mission characteristics of underground low-voltage networks - introducing techno-pedagogical (TP) method. International Jour- nal of Electrical Power and Energy Systems, 55, 253–260. Lucas, R., & Talukdar, S. (1978). Advances in finite element techniques for calculating cable resistances and inductances. IEEE Transac- tions on Power Apparatus and Systems, 3, 875–883. Marshall, J. S., Hines, P. D., Zhang, J. D., Minervini, F., & Rinjitham, S. (2013). Modeling the impact of electric vehicle charging on heat transfer around underground cables. Electric Power Systems Research, 97, 76–83. Orton, H. E.,&Samm,R. (1997).Worldwide underground transmission cable practices. IEEETransactions onPowerDelivery, 12(2), 533– 541. Papadopoulos, T. A., Chrysochos, A. I., & Papagiannis, G. K. (2013). Analytical study of the frequency-dependent earth conduction effects on underground power cables. IET Generation Distribu- tion and Transmission, 7(3), 276–287. Salles, M. B. C., Costa, M. C., Filho, M. L. P., Cardoso, J. R., &Marzo, G. R. (2010). Electromagnetic analysis of submarine umbilical cables with complex configurations. IEEE Transactions on Mag- netics, 46, 3317–3320. 123 Modeling of Power Cables with Arbitrary Cross Section: From Parameter Calculation to Electromagnetic Transients Simulation Abstract 1 Introduction 2 Parameter Calculation of Sector-Shaped Cables using FEM 3 Multiconductor Representation Using Modal Decoupling 4 Fitting a Modal Transformation Matrix for Cable Modeling in the Time Domain 5 Frequency-Dependent Modeling in the Time Domain 5.1 Frequency-Dependent Bergeron Circuit 5.2 Representation by State Equations 5.3 Cascade of Frequency-Dependent Bergeron Circuits 6 Time-Domain Simulations Using the Proposed Model 7 Conclusions Acknowledgements Appendix: Impedance and Admittance Parameters References