Investigations on the Developed Full Frequency-Dependent Cable Model for Calculations of Fast Transients †

A single-phase cable model based on lumped-parameters for transient in the time Abstract: The knowledge about the behavior of cables is substantial in cases of transients or in cases of faults. However, there are only a few models that are tailored to the current requirements for calculations of transient phenomena in three-phase cable systems. These models are based on complex structures. PI-section cable models with simple structures were previously qualiﬁed only for calculations in the frequency domain. A new full frequency-dependent cable model to simulate transient phenomena is introduced and validated. The model is based on lumped parameters with cascaded frequency-dependent PI-sections. For the implementation and the integration in simulation tools, it is important to investigate the impact of the PI-section parameters to the accuracy, the stability and the mathematical robustness. In this work, the impact of the frequency dependence of cable parameters, the length distribution and the number of PI-sections on the results of the developed three-phase cable model have been discussed. For simulations in the time domain, two algorithms have been presented to optimize the number of PI-sections based on a speciﬁed accuracy.


Introduction
Underground cables are an important part of electrical power systems in the HV, the EHV and the UHV level. Their application and the impact has increased with the growth of wind and solar energy to the total electricity generation [1,2]. The knowledge about the behavior in cases of transients is required for the correct design of cable lines and the overvoltage protection. The major difficulty for the modeling of cables is the complexity of the geometrical configuration and the electromagnetic coupling between phase conductors, cable shield and the earth. For this case, applicable models are required for the calculation of the electrical stresses. Effective models exist for load flow calculations and for calculations in the frequency domain. Older models were derived from approaches for overhead line calculations [3][4][5][6][7][8][9][10]. The results of these models have problems with result accuracy and the robustness. The incorrect modeling of geometrical condition and the incomplete interpretation of electromagnetic coupling are the reasons for the problems. To achieve a proper accuracy, a model with exact frequency-dependent cable parameters and with a good reproduction of wave propagation behavior is required.
In [11] a new full frequency-dependent three-phase PI-section cable model (3P-PI-Model) is introduced and validated. The frequency-dependent impedances in the model are determined with an improved sub-conductor method considering the skin and the proximity effect in the core, the shield and the earth [12,13]. The wave propagation in the model is based on PI-sections. The frequency dependence of the cable parameters is reproduced by means of vector fitting algorithm [14,15] and implemented through an RLC-network in a defined number of PI-sections. To decouple the three-phase cable system, a modal transformation is used [16][17][18].
For the implementation in calculation software, it is important to understand the behavior and the impact of the model parameter to the results. The significant problems of using PI-sections are their main resonance and the computational effort, which results from working with a large number of PI-sections. These problems can be influenced in the developed 3P-PI-Model through different parameters. The most important parameters, which have a direct impact, are the frequency dependence of cable parameters, length distribution of PI-sections and the number of PI-sections. To reduce the computational effort in the new model, two algorithms have been introduced for optimizing the number of PI-sections.

Wave Propagation Cable Model
The wave propagation in a multi-conductor cable or transmission line can be represented in the frequency domain with second-order differential Equations [19]: The longitudinal impedances Z and the shunt admittance Y L in the phase domain are calculated using the developed algorithms in [12,13]. The index "L" stands for lateral. In order to solve the differential Equations (1) and (2), a modal transformation technique has been proposed in [19]. As a result, the coupling between the cable conductors disappears and the cable system can be considered as decoupled single conductors. In the modal domain, Equations (1) and (2) can be rewritten as: where λ is a diagonal matrix which represents the eigenvalues matrix of Z Y L and Y L Z . For the transformation into modal domain the modal transformation matrices T i and T v are used. The relation between the matrices is given in [20]: Therefore, it is sufficient to calculate only one of them. In this paper, the results for T i will be calculated as in [17]. In order to perform the calculation in the modal domain, the quantities Z , Y L , V and I have to be transformed into this domain. From Equations (3) and (4), it can be derived: With simple matrix perturbation, the above equations are then rewritten as: By rearranging the terms and by using the Equation (5) in (8) and (9), the longitudinal impedance and the shunt admittance in the modal domain can be defined as: Both matrices in Equations (10) and (11) are diagonal. The frequency dependence of the elements of Y L,m cannot be neglected (as assumed for the elements of Y L ). This is related to the frequency dependent elements of T i . The elements of Z m are also strong frequency dependent. In order to take the frequency dependency of Z m and Y L,m into account, their diagonal elements will be approximated by mathematical functions using the vector fitting algorithm (VF) [14,15]. For example, the first diagonal element in Z m can be approximated using VF as: The accuracy of Z m fit depends on the number of the real terms N and the number of the complex terms M. Z m fit can be reproduced with an electrical network. In [21], a RL-network has been used. Using such a network allows to employ only the real poles and residues in the approximation (the complex terms are discarded). Although the fitting is achieved with a low order of approximation, the accuracy is affected. In this paper, an RLC-network shown in Figure 1 is introduced to regard all the poles and residues in the approximation including the complex terms. As can be seen, the network has been divided into two sections. The elements of the section I can be evaluated in a similar way as in [21]. The elements of section II are obtained as follows: The first diagonal element of Y L,m is similarly approximated. Furthermore, R 0 and L 0 are assumed to be zeros in Y L,m fit . After calculating Z m fit and Y L,m fit , they will be combined to form a developed PI-section. Figure 2 shows a representation for the developed PI-section after regarding the cable length, where v m,s and v m,r are the voltages in the modal domain at sending and receiving end respectively. One single conductor in the developed cable model can then be represented by cascading a number of PI-sections. Basically, the number of PI-sections depends on the cable length.  For the remaining diagonal elements of Z m and Y L,m , the same procedure is used to form the other conductors in the developed cable model. From the Equations (3)-(5), the voltages and the currents are transformed into the modal domain as follows: The Equations (14) and (15) are still in the frequency domain. However, since the simulations in the power system are necessary in the time domain, these equations have to be transformed into this domain. For the voltages in Equation (14), the convolution integral is applied as: Evaluation of this integral point by point is time consuming. By approximating the elements of T tr i using rational functions, an efficient recursive convolution method is used [7] and the above integral can be rewritten as a sum of exponentials: where b and l are the approximations parameters of the rational functions. The wave travel time τ for the impulse from one end of the cable to the other one is related to the imaginary term of the propagation function γ: In this relation, a 1 , a 2 and a 3 are constants that mainly depend on b, l and ∆ t [19]. The same procedure is used to transform the currents in Equation (15) into the time domain. In order to execute the time domain simulations, the differential equations for single conductors in the developed cable model have to be formulated. As shown in [21], the differential equations can be represented using state space techniques. A numerical Backward Euler rule of integration can then be applied to find the solution for the related state space equations [22]. To solve the cable equations with the rest of the network, which is always defined in phase quantities, the calculated voltages and currents in the modal domain must be transformed back to phase quantities. Figure 3 shows an overview for the realization of the 3P-PI-Model. Starting from defining Z and Y L , the cable model part is executed only once to provide the PI-sections for the calculation part. In this part, the computation is evaluated every ∆ t. For the first one, the vector of the voltages v is calculated using the network model and the vector of the currents i. By applying the recursive convolution between v and the approximation of t * i ( t * i ), the vector of the voltages v m is calculated. Using the PI-sections and v m , the state space equations are formulated and solved by means of numerical Backward Euler rule of integration to yield the current vector i m . The recursive convolution between i m and the approximation of t i ( t i ) provides the vector of currents i for the next ∆ t. Figure 4 shows a traditional PI-section without considering the frequency dependence of cable parameters for PI-section length of . As it can be seen, the main elements are divided into longitudinal (R and L) and shunt (C and G) parameters. In the study, constant permittivities are assumed for the insulation and semiconducting layers. It makes the shunt capacitance also frequency independent. The frequency dependence of the conductance is negligible [23,24]. As shown in [12,13], the strong frequency dependence of the cable longitudinal parameters makes it necessary to regard this dependency to simulate transient phenomena.

Impact of the Frequency Dependent Cable Parameters
Contrary to the traditional PI-Model, where the cable parameters are calculated at a certain frequency [25], the developed 3P-PI-Model includes the frequency dependence of the longitudinal impedance using an RLC-network (see Figure 1).
The importance of regarding the frequency dependence is well presented in Figure 5, where a comparison between the amplitude of the cable impedance at sending end with and without considering the frequency dependence of the longitudinal parameters is shown. The impedance response for both cases is calculated after connecting the receiving end of the cable to its characteristic impedance [25]. The cable type is given in Table 1. The calculations for the case without regarding the frequency dependence of the longitudinal parameters (traditional PI-Model) are carried out using the cable parameters at 50 Hz. For the case with taking the frequency dependence of the longitudinal parameters into account (developed 3P-PI-Model), the parameters are calculated using the proposed algorithm in [12,13]. As it can be seen, the amplitude of both models matches only at 50 Hz. The resonances of all PI-sections in both models accumulate and form a main resonance. This main resonance can be excited and therefore produce unwanted oscillations in the time domain, which do not exist in reality. This fact is a main problem when applying the PI-sections. However, as shown in this figure, for the same number of the PI-sections, the main resonance of the developed model is milder than one in the traditional PI-section model. In this case, the amplitude of the main resonance is about 10 times lower. The damping at the frequency of the main resonance is higher in comparison with damping in 50 Hz [13]. Therefore, considering the frequency dependence leads to a lower amplitude of the main resonance and as a result, to smaller unwanted oscillations.
Another difference between the results shown in Figure 5 is the frequency of the main resonance. The higher frequency of the main resonance in the developed 3P-PI-Model can be explained with considering the frequency dependence and therefore having a lower inductance in comparison to the traditional PI-Model [13,[26][27][28]. However, this fact does not directly have any considerable role in the results in the time domain. The first helpful effect here is that the main resonance is shifted thereby to a higher frequency with the higher damping. Secondly, a resonance at higher frequency is less likely to be excited during the calculations when the resonance is farther away from the frequency range of the simulated transients.
In sum, considering the frequency dependence can be a proper remedy for problems associated with applying the PI-sections.

Impact of the Length Distribution of PI-Sections
The accumulation of the individual resonance frequency of PI-sections leads to a main resonance (see Figure 5). This main resonance can be excited in the time domain, which will lead, depending on its amplitude, to unwanted oscillations. These oscillations can be influenced by distributing the length of PI-sections differently, which results in distribution profiles over the cable length .
The distribution profiles that give the best results are shown in Figure 6 for the cable length = 1 km and the number of PI-sections n PI = 20. The length of one PI-section is PI and its number is o PI (o PI = 1 ... n PI ).
In the case of uniform distribution, all PI-sections have the same length. In the parabola distribution, the length of each PI-section from the beginning to the middle of the cable is reduced. The section from the middle to the end of the cable is mirrored from the beginning to the middle of the cable. In the trapezoidal distribution, the first and last two PI-sections are shorter than the remaining PI-sections, which are evenly distributed.
The uniform distribution shown in Figure 6 is mainly used in the traditional PI-Model. Therefore, the impedance at the beginning of the cable when using the uniform distribution will be used as reference impedance in this section. Figure 7 presents the amplitude of cable impedance for the distribution profiles shown in Figure 6. The calculations are carried out for the 3P-PI-Model after connecting the receiving end of the cable to its characteristic impedance. The cable length is = 1 km. The number of PI-sections is n PI = 20.
As it can be seen, the distribution profiles shown in Figure 6 lead to the same impedance up to a certain frequency. The distinguishing characteristics of the impedances only relate to the amplitude and the location of the main resonance. The use of the parabola distribution results in a smaller amplitude of the main resonance compared to the uniform distribution. However, the main resonance of the parabola distribution has a lower frequency. The amplitude of the main resonance of the trapezoidal distribution is greater than that of the uniform distribution. However, the location of its main resonance has the highest frequency. This comparison leads to the conclusion that the distribution profile of PI-sections have an impact on the amplitude and the frequency of their main resonance.
As a compromise between the high frequency and the low amplitude of the main resonance of PI-sections, the uniform distribution is used for subsequent calculations.

Impact of the Number of PI-Sections
For investigating the impact of the number of PI-sections on the 3P-PI-Model's behavior, the distributed-parameters model, considering the frequency dependence of the cable parameters, is used as a reference. It should be noted that the reference model in the frequency domain can be developed without any simplifications. However, it is not directly applicable for the time domain calculations. To apply the reference model in these calculations, some simplifications are required [29][30][31].
For a single phase cable system, which is modeled through distributed-parameters, the current and voltage equations at the sending end are described as follows [19]: where γ = (R + jωL )(G + jωC ) According to Equations (20) and (21), if the receiving end of the cable is connected to the characteristic impedance, the cable impedance at the sending end will be the same as the characteristic impedance Z s = Z c (24) Figure 8 shows the amplitude of the cable characteristic impedance for the distributed-parameters model and the developed 3P-PI-Model for different numbers of PI-sections.
As shown, for a specific cable length and depending on the number of PI-sections, the behavior of the developed 3P-PI-Model will agree with the distributed-parameters model over a certain frequency range. By increasing the number of PI-sections, this range will be expanded.
For simulation of transients in the time domain, a large number of PI-sections has to be used in order to get a good accuracy. However, increasing the number of PI-sections leads to higher computational effort. This is caused by the higher dimension of the differential equations. Therefore, an optimization of the number of PI-sections is needed.

Analysis in the Time Domain
For modeling of cables with PI-sections, a minimum number of these sections has to be used. The number of PI-section depends mainly on the cable length and the simulated event. The longer the cable is, the more PI-sections are needed. In addition, more PI-sections are needed to simulate transient events (e.g., lightning strikes or switching operations) compared to steady state operations (e.g., during normal operation or failure).

Optimization of the Number of PI-Sections-First Algorithm
For optimizing the number of PI-sections in the developed 3P-PI-Model, a DC-Test on a three-phase cable with a length of = 2.5 km in trefoil formation has been used. The test is performed by charging the cable conductors up to 1 kV. After that, the charged conductors are connected to earth through the circuit breaker (s2) and the discharge current i core(A) at the sending end is simulated. The test configuration is shown in Figure 9. The cable parameters are given in Table 1.  Figure 9. Configuration for the DC-Test. Figure 10 shows the discharge current i core(A) with different number of PI-sections. As it can be seen, the slope of the current increases and the time delay of the current τ decreases. This is mainly due to the better replication of the higher frequencies in the current by increasing the number of PI-sections. However, the change in time delay dτ becomes smaller as the number of PI-sections increases. Above a certain number of PI-sections, this change can be neglected. Therefore, it is useful to set an upper limit on the number of PI-sections to optimize the computational effort. The upper limit of the number of PI-sections is determined using the algorithm shown in Figure 11. For a given cable length , the number of PI-sections n PI increases until the relative change in the time delay dτ r is less than a defined limit. The limit is initially set to 10 −4 . If the time delay is lower than this limit, then the loop is interrupted and the reached number of PI-sections is stored as the upper limit n PI L .
In this context, it should be noted that the overshoot of the current in Figure 10 (e.g., in the range 0 ≤ t ≤ 5 µs) is minimized and thus the amplitude of the current is also optimized by optimizing the time delay of the current. The reason is the dependence of the amplitude of the current on the characteristic impedance of the cable. This characteristic impedance is mainly calculated from the inductance and the capacitance of the cable, which in turn are included in the calculation of the time delay (see Equation (23)). Figure 12 shows the upper limit of the number of PI-sections depending on the cable length in the range 0.1 ≤ ≤ 10 km. For accurate simulation results, the number of PI-sections must be increased until the frequency of their main resonance is greater than the highest frequency of the event to be simulated. By increasing the cable length, the frequency of the main resonance of PI-sections will decrease [25]. In order to get the same accuracy of the calculations, the number of PI-sections has to be increased as it can be seen in Figure 12. For example, to simulate of transients by a cable length of 2.5 km the number of PI-sections is 45 and by a cable length of 10 km this number will increase to 52.
The proposed algorithm shown in Figure 11 defined the number of PI-sections in order to reach very good accuracy in the developed 3P-PI-Model, which results in using a large number of PI-sections. However, in many cases of simulation of transients a proper accuracy can be accepted. Therefore, the large number of PI-sections can be reduced and thus so can the computational effort.

Optimization of the Number of PI-Sections-Second Algorithm
The main idea of the second algorithm is to reduce the computational effort of the calculations in the developed 3P-PI-Model at the expense of the accuracy. Figure 13 illustrates the methodology of the second algorithm for the example shown in Figure 9. In this algorithm, the number of PI-sections is reduced to such an extent that a given error β in the time delay τ with respect to a reference value is not exceeded. The number of PI-sections will store under the reduced number n PI R . The reference value of the time delay is defined at upper limit by the number of PI-section n Π L , whereat the error in the time delay β is equal to 0 µs. By using the upper limit of the number of PI-sections shown in Figure 12 in the second algorithm the reduced number of PI-sections can be calculated for different errors in the time delay. Figure 14 shows the results for two specified errors. As it can be seen, specifying an error in the time delay leads to a considerable reduction in the number of PI-sections compared to their upper limit. For example, by a cable length of 2.5 km the number of PI-sections decreases from 41 to 12 by β = 1 µs and to 4 by β = 5 µs.
The impact of reducing the number of PI-sections on the simulation results will be illustrated using measured data. The measurements are performed using the same principle already described in Section 4.1. Figure 15 shows the measured current and the simulations with the developed 3P-PI-Model using the previously defined number of PI-sections. As it can be seen, using the upper limit of PI-sections (n PI = 41, β = 0) in the developed 3P-PI-Model results in a very good match between the measured current and the simulated one. This means that the wave travel time and attenuation are replicated correctly in the developed cable model. The match with the measured current is reduced when specifying an error in the time delay, whereby a worse match will occur by specifying a large error (β = 5 µs).
Reducing the number of PI-sections by specifying an error in the time delay affects the accuracy of the calculations. On the other hand, it results in a significant reduction of the computational effort. In this particular example, the dimension of the equations in the developed 3P-PI-Model are reduced from 6642 by n PI = 41 to 1944 by n PI = 12 and to 648 by n PI = 4. In addition, the computation time (MATLAB ® , typical PC-configuration) is reduced from 5 s by n PI = 41 to 2.5 s by n PI = 12 and to 2 s by n PI = 4.

Conclusions
In this work, detailed investigations of the behavior of a previously developed 3P-PI-Model are introduced. It has been proven that the frequency dependence of cable parameters, the length distribution of PI-sections and the number of PI-sections impact the frequency and the amplitude of the main resonance of PI-sections.
In addition, the number of PI-sections has a big impact on the computational effort in the developed 3P-PI-Model. Using a large number of PI-sections increases the accuracy of the simulation results. On the other hand, the large number of PI-sections leads to a high computational effort.
Two algorithms to optimize the number of PI-sections are developed. In the first algorithm, an upper limit of the number of PI-sections is defined in order to limit the computational effort in the developed 3P-PI-Model.
In the second algorithm, the computational effort of the calculations in the developed 3P-PI-Model is reduced at the expense of the accuracy. To achieve this, the number of PI-sections is defined (starting from the upper limit of PI-sections) after specifying an error in the time delay.
A comparison of the outcomes of the developed 3P-PI-Model using the upper limit of PI-sections with measured data has shown a good match, which validates the approach used in the calculating of the number of PI-sections. In addition, a proper accuracy of the simulation results with a large reduction in the computational effort can be reached by specifying a small error in the time delay.
In this work, the proposed first and second algorithms are used to define the number of PI-sections for simulating of fast transients. However, both algorithms can also be used by different operating states of the power system (e.g., switching operations).

Conflicts of Interest:
The authors declare no conflict of interest.