Novelty Underground Cable Model for Power System Transient Simulation

This paper presents an extension of the A-line model (Garcia-Sanchez et al., 2016) for underground cables that takes into account the frequency dependence and includes the conductance matrix; the new cable model is named the U-Cable model. This paper also presents the mathematical development that supports the new U-Cable model; thus, following the methodology, it can be observed how the cable parameters are related to the construction of the final mathematical model. This model has the characteristic that the fitting is generated by real poles in a natural manner without forcing vector fitting to change to only real poles, which makes the developed model extremely stable from a numerical point of view. The numerical Laplace transform (NLT) and PSCAD software are used for comparison purposes. The reason for the use of the NLT is that, if it is properly programmed, it always yields a reliable solution. For the comparisons, a short-circuited line is simulated first, followed by a matched line and then an open line with a unit step. A second simulation is performed by changing the unit step to a three-phase sinusoidal source. A comparison of the results obtained from the proposed model with those of the NLT and PSCAD provides support and reliability to this model because the differences observed, in most cases, are barely perceptible. Thus, it can be concluded that the presented model is an additional option to the existing models.


I. INTRODUCTION
The transmission of electrical energy by means of underground cable systems is used mostly in cities, principally where the risks of faults are important if aerial lines are used or when visual aspects are considered. Transient phenomena in cables are basically due to switching maneuvers and isolation failures. To reduce the risk of failure through enhanced cable design, adequate cable models must be used to simulate possible transient phenomena and their consequences. The simulation of transients in underground cables with frequency-dependent electrical parameters is performed with well-known different models, such as those proposed by Wedepohl and Wilcox [1], L. Marti [2], and the universal The associate editor coordinating the review of this manuscript and approving it for publication was Kai Sun . line model (ULM) [3], with the latter being the most accurate model.
In recent years, further investigations have been developed; for example, the frequency dependence of shunt admittances was treated in [4]. In 2009, Silva et al. [5] presented a new mesh domain model to simulate transients in underground cables. In 2019, Kocar and Mahseredjian [6] presented a new frequency-dependent model that belongs to the traveling wave class, which eliminates the numerical instability present in the ULM. A lossless line model was first developed for transmission lines based on the characteristics of partial differential equations for transient analysis of frequency-independent transmission lines and is presented in [7]. The electrical frequency-dependent parameters in the model were obtained by Ramirez et al. in 2001 [8]. Subsequently, the analysis of transmission lines with a method of characteristics VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ considering only the two transmission line end points was achieved by Escamilla et al. in 2013 [9]. The method of characteristics proposed previously takes as the starting point the transmission line equations developed by Radulet et al. [10]. This method requires transformation matrices to change from loop to phase quantities. A model based on the method of characteristics that does not require solutions at transmission line interior points and takes as the starting point the transmission line equations in the frequency domain is presented in [11]. Sánchez-Alegría et al. developed an extension of the model for field-excited transmission lines in 2019 [12].
The proposed model in this work, referred to as the U-Cable model, is a dual Norton circuit based on the method of characteristics presented in [11], which takes the transmission line equations in the frequency domain as the starting point and includes the frequency dependence of the conductance matrix. Moreover, the proposed method works directly in loop quantities, and subsequently, the solution in phase quantities is obtained via transformation matrices included in the proposed circuit. With this assumption, modal analysis is avoided in the solution procedure to solve the phase voltages and currents.

II. A-LINE MODEL INCLUDING CONDUCTANCE
Electromagnetic transient analysis in aerial transmission lines usually neglects the transversal conductance matrix, which in most cases is a valid assumption. However, when it is necessary to consider the conductance as constant, the method proposed in [11] can be extended.

A. TRANSMISSION LINE EQUATIONS
The transmission line equations in the frequency domain, considering the transversal conductance G independent of frequency, are defined as follows: where L G is the geometric inductance matrix, Z E (s) and Z C (s) are the earth return impedance and conductor internal impedance matrices, respectively, C G is the capacitance matrix, G is the conductance matrix, V(x, s) is the voltage vector, I(x, s) is the current vector, and s is the Laplace variable. The earth return impedance and the conductor's internal impedance can be grouped in a penetration impedance Z , as defined in [11].
Transforming (1) to the time domain: where ⊗ denotes a convolution operation, which can be computed by substituting Z with a rational function with poles, residues, and a constant matrix. Solving the convolution using the auxiliary differential equation method, we obtain, for the case of (2a): where p n and k n are the poles and the corresponding residue matrix, respectively, and r k is a real constant matrix. Moreover, from (3b): Expressions (2b) and (3a) are the time domain transmission line equations that should be solved. Equation (3a) is an inhomogeneous partial differential equation, where the excitation function (3b) at time t is given in terms of the delayed current values (4).

B. METHOD OF CHARACTERISTICS
The characteristic impedance of an ideal multiconductor transmission line is defined as follows: Applying the method of characteristics as proposed in [11] and restricting the solution along the curves shown in Fig. 1 and defined by: (2b) and (3a) can be written as dv dx Solving (7a) using the backward Euler method between transmission line ends A at t-τ and B at t yields: where and U is the identity matrix. Solving (7b) using the backward Euler method between the transmission line ends A at t and B at t-τ yields: Considering the currents entering the line at both ends, and substituting (11) into (8), we obtain, in a compact form, The vector h B represents a history current source connected to the line-receiving end.
Substituting (11) into (10), we obtain, in compact form, The vector h A represents a history current source connected to the line-sending end. Equations (12) and (13) represent the dual Norton model of aerial multiconductor transmission lines, taking into account the conductance matrix. The equivalent circuit is shown in Fig. 2. If the conductance matrix is neglected, these equations turn into those proposed in [11].

III. MODEL FOR UNDERGROUND CABLE
In this section, the method of characteristics previously exposed is extended for application to underground cables considering the frequency dependence of the conductance matrix. This model is called the U-Cable model.

A. CABLE LINE EQUATIONS
The transmission line equations for underground cables in the frequency domain, splitting the geometric part (as shown in Appendix A), are defined as follows: where the subscript l indicates that the equations are in loop quantities. The conductance matrix G with frequency dependence is included. These equations in the time domain are expressed as follows: Solving the convolutions was performed using the procedure described in [11], yielding: and from the curve fitting technique, the subsequent states are obtained, where the pairs (p n , k n ) are the poles and the corresponding matrix of residues, r k is a real constant matrix, and N is the order of the rational fitting of z . In contrast, the pairs (q n , m n ) are the poles and the corresponding matrix of residues, G k is a real constant matrix, and N is the order of the rational fitting of g . Furthermore, The inhomogeneous partial differential Equations (16) are the telegraph equations that should be solved for the case of underground cables. The excitation functions E and F at time t of these equations are given in terms of the delayed values of the currents and voltages, respectively, in (19).
As seen in (15), the model proposed in this work requires fitting of the penetration impedance Z and conductance G to solve the convolutions, which is achieved using vector fitting.

B. CHARACTERISTICS METHOD FOR CABLES
To apply the method of characteristics proposed in [11], it is required that the product of L G and C G yields a diagonal matrix so that working with a decoupled system is possible. As this is accomplished in a straightforward manner when working with the geometric inductance and capacitance parts in loop quantities, for an ideal cable, the characteristic impedance is defined as follows: and subtracting (21) from (16a) gives so, it is easy to show that where U is the identity matrix. Substituting the definitions (23) into (22), Each element of λ in (23) is the inverse of the speed of the wavefront in the corresponding i-th dielectric material: Thus, for different dielectric materials, this is, Fig. 3.
Restricting the solution of (24) along the curves defined by (25), (24) can be written as:

C. ANALYSIS OF THE NUMERICAL IMPLEMENTATION
Solving (26a) using the backward Euler method between the transmission line ends A at t − τ and B at t yields: solving from (27) for the current at B, Applying the same integration method to (26b) between the transmission line ends A at t and B at t − τ we obtain: solving for the current at A, Considering that the positive direction of the currents enters the line at its ends, as in (11), (30) becomes: Or in compact form: The vector h A is given in terms of past values of loop voltages and currents; therefore, it represents a history current source connected at the sending end of the cable. Now using (11) in (28), gives: Or in compact form: is a vector h B that represents a history current source connected to the receiving end of the cable. Equations (31) and (32) represent the dual Norton model of the underground cable in loop quantities, which includes the frequency dependence of the conductance matrix.
These equations consider the general behavior of the underground cables and transmission lines. Considering the case where the conductance matrix does not depend on the frequency, the term F in Equations (31) and (32) becomes zero, which turns the equations into (12) and (13). Additionally, the term η 1 becomes the unity matrix for the case neglecting the conductance matrix, as mentioned above; thus, the equations proposed in [11] are obtained.
The equations are expressed in loop quantities. To work properly with the proposed model in its Norton structure connected to external elements so that core, sheath and armor voltages and currents can be used with the ground as a reference, transformations to the phase domain are necessary. Nevertheless, applying the transformations shown in Appendix A turns the product of L G and C G into a full matrix, which makes it necessary to use additional transformations to obtain a decoupled system so that the proposed model can be applied.
For converting phase to loop quantities, the following transformation matrices are required [13]: Thus, substituting (33a) into (31) and (32), we obtain: solving for the currents, where the history currents, now in phase quantities, are The loop currents I A l,t and I B l,t can be represented as currentcontrolled sources via the transformation matrix T l in (33). For each i-th loop current where N is the number of loop currents, and the rows and columns of matrix T i are indicated by the subscripts i and k, respectively. Moreover, for the phase current vector, the subscript k represents the k-th row element.
In the case of voltages, from (34), The voltages v A φ,t and v B φ,t can be represented as voltagecontrolled sources via the transformation matrix T −1 v ; for each i-th row, we have: that can be expressed as: where N is the number of conductors and the rows and columns of matrix T v are indicated by the subscripts i and k, respectively. Moreover, for the phase voltage vector, the subscript k represents the k-th row element, and The corresponding circuit with voltage-and currentdependent sources is shown in Fig. 4 for a cable with a core, sheath, and armor. Note that the central part of the circuit corresponds to loop voltages and currents, where the proposed method (U-Cable model) in its Norton representation is applied. In contrast, the outer part of the circuit corresponds to the phase quantities. The core, sheath, and armor circuits are interconnected by means of dependent voltage and current sources.
In the circuit, h A and h B are the current vectors that contain the core, sheath, and armor delayed currents. Furthermore, in the admittance matrix (G L η 1 ) i , the subscript i corresponds to the i-th row of the matrix.
Subscripts 1 and 2 in the resistances are used to indicate that the corresponding element is connected at the sending and receiving ends, respectively. For multiple cables, the coupling between external conductors of each cable through the ground must be considered. This coupling, owing to the outer loop currents (Fig. 5), is considered in the following equations for a 3-cable system with core (c), sheath (s), and armor (a), where (x, t) has been omitted for the sake of simplicity: The numeric subscripts indicate the number of cables, and the subscript ''a'' designates the armor voltages and currents. The convolution procedure of the mutual earth impedance z gm,i−j with the armor currents i j,a is accomplished using the auxiliary differential equation method as in [11]. This procedure for each mutual coupling of the outer currents between cables i and j yields: where the pairs (p gm , k gm ) are the poles and the corresponding residues, respectively, r gm is a real constant matrix, and N is the order of the rational fitting of z gm .

IV. APPLICATION EXAMPLES
The proposed model was proven using two different application examples. The first consists of three-phase underground cables with a unit voltage step applied; in this example, a comparison is made against the results from the NLT and PSCAD. The second is a three-phase underground cable system with a three-phase source.

A. THREE-PHASE CABLE WITH UNIT CURRENT STEP
The analysis of the proposed model in three-phase cables was carried out with a unit current step injection at the sending end of each core. The cable configuration is shown in Figs. 5 and 6; the corresponding data are shown in Tables 1 and 2. The sending and receiving ends of the sheath and armor are connected to the ground by means of resistances of 1 × 10 −6 . The conductance matrix is considered independent of the frequency. First, in the receiving end, the core was left short-circuited with a 0.001 resistance. As a second test, the receiving end of the core was matched with the characteristic impedance, which was 26.585 . Finally, the receiving end of the core was open with a resistance of 1000 . The voltages obtained with the proposed method were compared with those obtained from the NLT and PSCAD. Figs. 7, 8 and 9 show the voltages at the sending end for the three conditions. Figs. 10, 11 and 12 show the voltages at the receiving end.   The analysis of the graphs shows that, for the case of a unit step source, the results obtained with the U-Cable model and the NLT are similar, and a small difference from PSCAD is present. Only the results of one of the phases of the system are presented; owing to the practically identical results of the other phases, their results are omitted.

B. THREE-PHASE CABLE WITH AC SOURCE
The following simulations consist of using a sinusoidal source instead of a unit step source; the line conditions are maintained. In the first case, a short-circuited line is simulated at the receiving end. Fig. 13 shows the voltage in the core for the three cables, and Figs. 14 and 15 show the voltages in the corresponding sheaths and armors, respectively, all of them at the sending end.        When analyzing the results obtained, not only in the core but also in the sheath and armor of each cable, it is noted that the source has a phase to neutral voltage of 132.79 kV   or 230 kV between lineswith an internal source resistance of 1.587 . In addition, it is noted that there is voltage damping, so an 80 V phase to neutral is obtained, which is due to the short-circuited end with a resistance of 0.001 . Even though a more realistic value of resistance could be used,  Connecting the receiving end with its characteristic impedance, which in this case corresponds to 26.585 , gives the next test performance. Figs. 19 and 20 show the voltages at the sending and receiving ends, respectively. Analyzing the results, it is noticed that the waveforms have practically no reflections. The voltages due to mutual couplings are quite similar to those obtained in previous simulations; therefore, they have been omitted because they do not yield additional information.
Finally, the line is analyzed as an open circuit with a resistance of 1000 connected at the receiving end; thus, compared with the characteristic impedance, it behaves as an open line. The behavior at both ends is shown in Figs. 21 and 22. A higher value for the resistance can be used, but the results do not change at all except for the fact that if the resistance is higher, or in other words, there is a greater difference between the characteristic impedance and the connected resistance, the error in the numerical convergence is larger. Thus, there is no reason to justify the use of higher resistance.

V. DISCUSSION
Electric power network analysis and simulation in the transient state is necessary to perform, including the frequency dependence of all its elements; the traditional method of simulation uses the numerical Laplace transform (NLT), in which a sequential simulation is not possible because one of the relevant parameters necessary for its implementation is the observation time. For this reason, different techniques in the time domain have been created; the universal line model (ULM) is the most well-known. However, the search for alternatives is open because it does not seem feasible that a specific numerical technique meets all the requirements and fits for all the test cases and in all senses. There is an effort to develop a suitable model for each case and need using the FdLine model and techniques such as finite element modeling, neural networks and finite differences, among others. In this sense, the proposed U-Cable model was developed for practical purposes, including comparisons with solutions obtained from NLT and PSCAD. Naturally, this model has its advantages and disadvantages without detracting from any of the previous techniques.
Although it is known that no numerical method can be a reference for another, which is the main reason for not including the results of the other existing methods in this work, it is also true that there are techniques that, if they are properly implemented, can always obtain good results. When there are no restrictions regarding the execution time, the results obtained are accurate, as in the case of the NLT. That is, the natural technique to simulate a system with frequency dependence is the NLT. Unfortunately, this technique has limitations that do not allow the simulation of relatively small systems. For example, if there is a three-phase network with 4 generators, 4 loads, and 10 power lines, then the equation system is 54 × 54, which has to be inverted for each frequency with an accumulation of error. Because each of the partial results participates in each point in the time domain, the local errors of each inverse when passing to the time domain become a global error.
This would be very easy to verify because we can purposely make an inverse with an error, for example, of 50% in the frequency domain and see how this error affects all points in the time domain. In conclusion, the NLT is the best reference point in the case of simulations with few nodes, but it is not the best tool for network simulations, which is a limitation that restrains the use of the NLT and forces the search for alternatives such as those previously mentioned, among others, and those that are explored day-by-day. Thus, the technique proposed in this work adds to the accumulation of models in the test and study phases.

VI. CONCLUSION
In this paper, a novel methodology for simulating underground cables with frequency dependence is presented. For the parameter fitting stage, a technique called vector fitting (VF) is used, which is well proven [14]- [16]. The results obtained were compared qualitatively and quantitatively with those obtained with the NLT and the PSCAD software.
The simulations are performed with two types of sources: DC and AC with the receiving ends in short circuit, matched, and open lines. The main findings are as follows: 1. -When the DC source is connected, the voltages at the sending and receiving ends of the core match those of the NLT and PSCAD software. At first glance, the maximum deviation can be observed in the short-circuit simulation.
2. -For the case of voltages in both the armor and the sheath, there are more notable differences, but it must be taken into account that the amplitudes are near 10 −7 pu. These differences can be attributed to both fittings and the way the model couples the cables.
3. -For the case of the AC source, we first present the results related to a short-circuit in the receiving ending, where it can be seen how the voltages at both ends are quite similar to those obtained with the NLT and the PSCAD software to the point that they are superimposed and not fully distinguishable. The time scales serve to show how in the sending end side there are a few oscillations and the voltage stabilizes at approximately 100 kV, and on the receiving ending side, where there is a sustained short-circuit, there is a voltage of 80 V. In the case of both sending and receiving couplings, there are similar amplitudes and shapes because they are coupled with the same resistance.
4. -The following simulation corresponds to one of the matched lines where it can be seen that there is a small reflection on the remote side because the line is initially de-energized but immediately the voltages stabilize. Only those voltages related to the core are presented because the sheath and armor voltages do not change in connection, so the results are practically identical to the previous ones. 5.
-For the open line condition, the connectivity of the sheath and armor is maintained, so the results continue to be similar to the previous ones, which is the reason why only those voltages on the core are presented. These results show how the reflections attenuate until they disappear, and the resulting voltage behaves as if the line is matched.
6. -From the analysis of the results, it can be concluded that the proposed methodology is adequate for the simulation of underground cables with frequency dependence in the time domain.
7. -The fitting stage was performed naturally using only real poles. The reason why the results are not presented is because they are not conclusive; that is, a study is required for all physically feasible cable configurations so that it can be determined that the model can be fitted with only real poles in all cases. However, these results are promising. Seven real poles were used in the test cases.

APPENDIX A
When working with underground cables (Fig. 23), the electrical parameters are obtained by analyzing the voltages and currents in loop quantities, as shown in Fig. 24. The telegraph's equations of voltage in the frequency domain, omitting (x, s) for simplicity, are [13]: where Z G is the geometric impedance, Z bb is the outer surface impedance, Z aa is the inner surface impedance, and Z ab is the transfer impedance for the core, sheath, armor, and earth, identified by the subscripts C, Sh, A, and E, respectively. These cable parameters are calculated using the formulae proposed by Wedepohl and Wilcox [1], whereas the earth parameters are obtained using the Saad-Gaba-Giroux approximation [17]. Note that in (A1), the geometric part is separated from the cable and ground impedances. On the other hand, the equations of current are:  The matrix form of the equation system in loop quantities is shown below: where L G and C G are diagonal matrices given by and For working with quantities referred to ground, that is, the core, sheath, and armor voltages and currents as reference, it is necessary to change from loop to phase quantities via the following definitions obtained from Fig. 7 [13]: Thus, the voltage equation system of (A3) in matrix form considering (A7) is: where where Z s1,1 = Z 1,1 + Z 2,2 + Z 3,3 + 2Z 1,2 + 2Z 2,3 (A11a) Z s2,2 = Z 2,2 + Z 3,3 + 2Z 2,3 (A11b) Z s3,3 = Z 3, 3 (A11c) Z s1,2 = Z s2,1 = Z 2,2 + Z 3,3 + Z 1,2 + 2Z 2,3 (A11d) Z s1,3 = Z s3,1 = Z s2,3 = Z s3,2 = Z 1,1 = Z 3,3 +Z 2,3 (A11e) In addition, the current equation system in (A2) in matrix form considering (A7) is: When conductance G is included in the analysis, the admittance matrix must be split as follows: (A14) Note that both the geometric impedance and admittance matrices are diagonal but the resulting product of the converted matrices to phase quantities is not. Hence, the proposed method cannot be implemented directly by transforming the telegrapher's equations into phase quantities without employing modal transformations. However, note that the product of the geometric matrices in (A4) yields a diagonal matrix such that the proposed method can be implemented with loop quantities.