Enhancing the frequency-domain calculation of transients in multiconductor power transmission lines

Electric Power Systems Research 122 (2015) 56–64

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

Enhancing the frequency-domain calculation of transients in multiconductor power transmission lines Andreas I. Chrysochos, Theofilos A. Papadopoulos, Grigoris K. Papagiannis ∗ Power Systems Laboratory, School of Electrical and Computer Engineering, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece

a r t i c l e

i n f o

Article history: Received 7 September 2014 Received in revised form 24 December 2014 Accepted 26 December 2014 Keywords: Electromagnetic transients Fast Fourier transform Frequency-domain modeling Levenberg–Marquardt method Numerical Laplace transform Shape-preserving piecewise interpolation

a b s t r a c t In this paper numerical improvements are proposed to enhance the accuracy and efficiency of frequencydomain transient transmission line modeling. The adopted procedure is based on a robust eigenvector tracking algorithm and an accurate shape-preserving interpolation routine, formulating the highly resonant nodal admittance matrix from a minimal number of frequency data. The proposed method is integrated with the numerical Laplace transform and its performance is demonstrated in a configuration, consisting of an underground cable and an overhead line. Results are validated with the corresponding obtained from a frequency-domain model with regular sampling, revealing the high accuracy of the proposed formulation. The computational efficiency of the proposed numerical improvements is highlighted by comparisons of the execution times in cases of high sampling rates. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The development of simulation models for the accurate calculation of electromagnetic transients in overhead lines and underground cables has been considerably increased again in recent years. A large number of simplified and advanced models have been proposed, which can be generally classified into two main categories according to the solution method adopted in the calculation routine. In time-domain modeling the solution is based on a simple rule of integration, mainly exploiting the benefits of the recursive convolution technique [1]. This is significant, since such models can be easily incorporated in EMTP-like programs. These models can be further divided into modal- [2,3] or phase-domain [4], depending on whether the transfer propagation and the characteristic impedance/admittance matrices are expressed and fitted in their diagonalized or full matrix form, respectively. On the contrary, in frequency-domain models, calculations take place solely in the frequency-domain, while the voltages and currents are converted back to the time–domain with the explicit use of the inverse Fourier or Laplace transform [5–7]. Although

∗ Corresponding author at: Power Systems Laboratory, School of Electrical and Computer Engineering, Aristotle University of Thessaloniki, P.O. Box 486, 54124 Thessaloniki, Greece. Tel.: +30 2310996388; fax: +30 2310996302. E-mail address: [email protected] (G.K. Papagiannis). http://dx.doi.org/10.1016/j.epsr.2014.12.024 0378-7796/© 2015 Elsevier B.V. All rights reserved.

frequency-domain models were initially used in linear networks, considerable efforts have been made to include non-linear and time-varying elements, network configuration modifications and non-equidistant sampling schemes [8–15]. This class of models is usually used as stand-alone transient analysis tools or as backup options for the validation of time-domain models in certain cases. One of the major drawbacks in frequency-domain modeling is the number of samples required to resolve the transient responses [13], which is explicitly determined by the frequency sampling interval. The latter can be extremely small in cases of transmission network simulations, where long observation times and extended frequency ranges are required due to different lengths and resonance frequencies [14]. As a result, an excessive number of samples is required leading to high computational costs in the preparatory calculation stage of the per-unit-length (PUL) parameters and propagation characteristics. The same issues can also occur in repetitive studies [16], such as statistical analysis, where multiple simulations of the examined system are performed by varying specific parameters of interest. Scope of this paper is to present numerical improvements for the accurate and efficient frequency-domain calculation of transients in multiconductor power transmission lines (TL). More specifically, a backward eigenvector tracking algorithm and a shape-preserving interpolation technique are proposed for the formulation of the highly resonant nodal admittance matrix from a relatively small number of initial frequency data, which is predefined by the user. The proposed numerical improvements, combined with the

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64

utilization of the numerical Laplace transform (NLT) [6], provide high accuracy in the simulation of any type of transients and great computational efficiency, even in cases requiring high sampling rates or when time-consuming and complex integral forms of earth representation must be used for the calculation of the PUL parameters. The proposed numerical improvements are applied to a power line configuration, consisting of an underground cable and an overhead line in cascade connection, highlighting the robustness of the eigenvector tracking algorithm and the high precision of the interpolation procedure. The accuracy of the resulting nodal admittance matrix is thoroughly investigated for variable number of initial frequency data. The obtained transient responses are further validated with the corresponding derived from a standard NLT model, where all frequency samples are directly calculated. Finally, the computational costs are compared, revealing the high efficiency of the proposed method, especially in cases where complex representations of the imperfect earth are assumed.

2.1.1. Step 1 Assuming the regular sampling scheme of [6] and P equally spaced samples in time- and frequency-domain, the standard method starts with the derivation of the required angular frequency and time intervals ω, t from the user-defined total observation time T and maximum examined angular frequency ˝. t = /˝

(1)

ω = 2/T = 2˝/P

(2)

2.1.2. Step 2 The N × N PUL series impedance and shunt admittance matrices in frequency–domain, Z and Y respectively, are calculated at each frequency defined by the interval ω. Modal-domain analysis for multiconductor lines is applied to the matrix product Y Z at each discrete frequency to calculate the modal propagation characteristics, namely the propagation constant matrix  as well T

as the corresponding eigenvector matrices Ti and T v = [T i −1 ] [17]. The complex and frequency-dependent eigenproblem can be solved using the simple and fast-converging QR diagonalization algorithm [18,19]. Any possible numerical eigenvector switchovers are treated in the final formulation of the nodal admittance matrix (Step 3), which is expressed in the phase-domain [20]. 2.1.3. Step 3 The power line can be represented by the exact N-phase PIequivalent circuit of Fig. 1 [15]. The admittance matrix Y PI−eq of (3), connecting the two nodes of Fig. 1 and also including the termination conditions YS , YR at both ends, is first defined. Y sr + Y sh + Y S

−Y sr

−Y sr

Y sr + Y sh + Y R

  =

A + YS

B

B

A + YR



(3)

where: A = Z

−1

B = −Z

· T v · diag{} · diag −1 {tanh( · )} · T Ti

 −1

⎡

Y 11

Y M1

In this section the standard NLT model is summarized [6], presenting the steps for the calculation of the time-domain transient responses on an N-conductor power line with length .



nodal analysis approach. This matrix is calculated at each discrete frequency of the total number P defined by the adopted sampling scheme in (1) and (2), and is of order L × L with L = N × M, where M is the number of nodes in the PI-equivalent circuit.

⎢ . ⎣ ..

2.1. Overview

Y PI−eq =

Fig. 1. N-phase PI-equivalent circuit terminated at both ends.

Y nd = ⎢

2. Standard NLT transient model

· T v · diag{} · diag −1 {sinh( · )} · T Ti

(4) (5)

Generalizing (3) for topologies with more than two nodes, the nodal admittance matrix Ynd of (6) is formulated by applying the

57

···

Y 1M

..

.. .

.

...

⎤ ⎥ ⎥ ⎦

(6)

Y MM

In (6) the symmetric Yii , of order N × N, represents the self and mutual coupling of the N conductors for node i over the examined frequency range, whereas the symmetric Yij , of the same order, describes the mutual coupling between nodes i and j for all N conductors [15]. 2.1.4. Step 4 Utilizing the algorithm of the NLT [6], the voltage and current sources are transformed in the frequency-domain with the same frequency interval used in the nodal formulation of Step 3. Therefore, the node voltages End and the injected node currents Ind , both of order L × 1, can be related with the generalized nodal admittance matrix Ynd , resulting in (7). I nd = Y nd · E nd

(7)

2.1.5. Step 5 Eq. (7) is reduced by eliminating the rows and columns corresponding to the known nodes. The remaining voltages can be easily calculated by (8) with Erd containing the unknown node voltages, whereas Ird and Yrd are known. Results are finally converted back to the time–domain, implementing the algorithm of inverse numerical Laplace transform (INLT) [6]. E rd = Y −1 · I rd rd

(8)

2.2. Sampling scheme considerations According to the regular sampling scheme adopted in the calculation routine of the standard NLT model, there are several options to determine T, ˝ and ω, t. In most cases, once T is selected, ω is defined by (2). Next, depending on the type of transient under analysis, ˝ is selected according to the required bandwidth that has to be considered. Thus, t and P are calculated by (1) and (2), respectively. From the considered sampling scheme it is evident that in the worst case, i.e. for long simulation times and for wide frequency ranges, a large number of equidistant samples is required [14]. As an example, 104 samples (without counting the conjugate symmetric part) are needed for a total simulation time of 10 ms and a frequency range between 0 Hz and 1 MHz, resulting in an equal number of calculations for the PUL parameters and propagation characteristics in Step 2. This increases significantly the computational burden, as it is the most time-consuming process [13], especially when imperfect

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64

υ

58

Fig. 2. Underground cable cross-section.

earth formulations, involving numerical calculations of complex semi-infinite integral forms, are considered [21]. 3. Proposed numerical improvements A solution to this issue can be the use of numerical interpolation for the calculation of the complex nodal admittance matrix from a minimal set of initial frequencies, defined as  = {f1 , f2 , . . ., fm } . In this paper, a shape-preserving piecewise interpolation technique [22,23] is applied to Z ,  and Ti that yield A, B (4) and (5), in order to calculate Ynd at each discrete frequency defined by the interval in Step 1. This approach is designed to respect the curvature of the interpolated data, minimizing the error compared to other simple polynomial fittings. The calculation of  and Ti is performed at the same set of frequencies , by utilizing an eigenvector tracking algorithm to avoid any abrupt discontinuities. The proposed method adopts the robust computational scheme of [24], resulting in smoother eigenvalues and eigenvectors without errors even for coarse sampling rates. As a result, an indirect but very efficient calculation of the complex admittance matrix is performed from a minimal set of frequencies  without loss of accuracy. On the other hand, a direct numerical interpolation of Ynd elements from the same set  would result in significant computational errors, due to the highly resonant behavior of Ynd in the frequency-domain. 3.1. Eigenvector tracking algorithm Due to the proposed interpolation procedure for the indirect calculation of Ynd the propagation constant matrix  and the eigenvector matrix Ti calculated in Step 2 must be smooth functions of frequency without switchovers. In this paper, a recently developed backward eigenvector tracking algorithm based on the Levenberg–Marquardt (LM) method is selected [24]. The algorithm solves an equivalent real-valued formulation of the complex eigenproblem, providing smooth eigenvalues and eigenvectors in all line and cable configurations, regardless of the frequency sampling rate and without any calculation errors [24]. The improved performance of the adopted algorithm over other implementations and especially the Newton–Raphson (NR) method [25] is demonstrated in the three-phase, medium-voltage, single-core cable system of Fig. 2, which is buried 1-m under the ground surface in a flat formation. Assuming the homogeneous earth formulation of [26], the tracking ability and robustness of the LM method is examined by varying the sampling rate. Instead of 104 samples that would be normally required by the standard NLT for the frequency range of 0–1 MHz and a total time of 10 ms, the six modal propagation velocities (ground mode #1,

Fig. 3. Modal propagation velocities of the underground cable calculated by the LM and the NR methods with 3 and 10 points per decade.

intersheath modes #2–#3, and coaxial modes #4–#6) are calculated in Fig. 3 for the LM and the NR methods, assuming a set  of only 3 and 10 points per logarithmic decade. Although the former sampling rate is seldom used in most practical applications, it is a good indicator of the numerical stability of each algorithm. The NR method with 3 points per decade fails to calculate the propagation velocities of both intersheath modes #2–#3 above certain frequencies, as highlighted in the figure. The LM method preserves its tracking ability, even in the case of this large frequency interval, due to the enhanced line-search criterion that is incorporated in the calculation routine [24]. The LM method provides also smoother eigenvectors compared to the NR method in cases of both high and low sampling rates [24], which is a rather significant attribute for the accurate numerical interpolation of the corresponding propagation characteristics. 3.2. Interpolation technique After the diagonalization procedure, Z ,  and Ti calculated at the minimal set  are interpolated using piecewise functions and are further evaluated at each discrete frequency defined by the interval of Step 1. This procedure is based on a modified version of the Fritsch-Carlson method [22,23] and is separately applied to the real and imaginary part of each matrix element. The resulting new variables are denoted with the superscript f, which express the extension and evaluation in the desired spectrum. Assuming the breakpoints (xi , yi ) for each matrix element at the set of frequencies , i = 1,2,. . .,m, the interpolated function L(x) at the i-th interval xi ≤ x ≤ xi+1 is expressed in terms of the local variables s = x − xi and hi = xi+1 − xi ,



L(x) =

3hi s2 − 2s3

+



h3i s2 (s − hi ) h2i



yi+1 +



h3i − 3hi s2 + 2s3

di+1 +

h3i s(s − hi ) h2i

2



yi

 di

(9)

satisfying the following four interpolation conditions: yi = L(xi ),

yi+1 = L(xi+1 )

(10)

di = L (xi ),

di+1 = L (xi+1 )

(11)

The yet unknown slopes of (11) are determined in a way that function L(x) does not locally overshoot the initial data. Therefore,

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64

59

(b)

(a)

Fig. 4. (a) Magnitude and (b) angle part of five elements of the PUL series impedance matrix for the underground cable.

the slope di for all interior points (i = 2, . . ., m − 1) is defined by the sign comparison of the differences ıi−1 , ıi , given in (12) ıi−1 =

yi − yi−1 , xi − xi−1

ıi =

yi+1 − yi xi+1 − xi

(12)

If ıi−1 and ıi have opposite signs or if either of them is zero, the slope di is set to zero since (xi , yi ) is a local extremum. On the other hand, in the case of the same sign, the slope di is calculated as the weighted harmonic mean of (13), with weights determined by the lengths of the two consecutive intervals in (14). w1 w2 w1 + w2 = + di ıi−1 ıi

(13)

w1 = 2hi + hi−1 ,

(14)

w2 = hi + 2hi−1

Finally, the first and last slopes at the endpoints (d1 , dm ) are calculated by the one-sided formulas of (15), respectively, checking also the signs of the corresponding differences in order to avoid abrupt overshoots [23] and possible dc bias errors at transient response tails. d1 =

(2h1 + h2 )ı1 − h1 ı2 , h1 + h2

dm =

(2hm−1 + hm−2 )ım−1 − hm−1 ım−2 hm−1 + hm−2

(15)

In Fig. 4 the magnitude and angle of five elements (i,j) of the f 6 × 6 PUL series impedance Z  are shown for the cable configuration of Fig. 2. The interpolated curves, obtained by the interpolation technique of (9)–(15) and plotted with frequency interval 1 Hz, are compared to the corresponding data, derived at set  with sampling rate 10 points per decade, showing high accuracy. Due to the cable cross-section symmetry, the rest elements of the 6 × 6 PUL impedance matrix follow the same pattern and thus are not illustrated. Furthermore, in Fig. 5a and b the propagation characteristics of  f are presented, as decomposed into the six propagation modes, namely ground mode #1, intersheath modes #2–#3, and coaxial modes #4–#6. In Fig. 5c and d the magnitude and angle part of the first column (i,1) of the 6 × 6 modal-domain transformation matrix f T i is plotted, corresponding to the ground eigenvector #1.

From the obtained results it is evident that the integrated interpolation preserves the shape of the data at set  and respects the corresponding monotonicity without overshoots and abrupt oscillations. This is mainly attributed to the particular selection of the slopes (12)–(15), which results in a continuous first but not second derivative and implies a discontinuous but also shapepreserving curvature. Therefore, a high degree of precision is generally achieved. Moreover, this technique is preferred over simple polynomial fittings, since the error is significantly smaller, avoiding also the problem of Runge’s phenomenon [28]. 3.3. Incorporation in frequency-domain transient model The inclusion of the proposed numerical improvements in the NLT transient model of Section II is shown in Fig. 6 by means of a workflow, highlighting the Modified Steps for the calculation of time-domain transient responses. The proposed methodology can be integrated with any discrete Fourier transform algorithm preserving its numerical advantages, since it actually enhances the numerical efficiency in the formulation of the nodal admittance matrix and does not modify the discrete algorithm used. Thus, the Fast Fourier Transform can also be adopted, although it may be less accurate compared to the NLT implementation. This is attributed to Gibbs oscillations (truncation error) and to aliasing phenomena (discretization error), caused by the absence of the window function and the damping factor, respectively [6,27]. Furthermore, the proposed numerical modifications are incorporated in the regular sampling scheme, which is rather simple and flexible, since it can be used in any type of transient without prior knowledge of their behavior. The proposed improvements can also be integrated with any other sampling scheme, such as the more efficient odd sampling, where Gibbs oscillations are even more attenuated due to the double upper frequency limit [6]. 4. Results 4.1. Numerical accuracy The combined performance of the proposed numerical improvements is first demonstrated in the calculation of the nodal

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64

α

υ

60

Fig. 5. Modal (a) attenuation and (b) propagation velocity. (c) Magnitude and (d) angle part of the corresponding ground eigenvector for the underground cable.

Fig. 7. Overhead line cross-section.

Fig. 6. Workflow of the NLT transient model with the proposed numerical modifications.

admittance matrix for the 1-km cable configuration of Fig. 2 and the 50-km three-phase, medium-voltage, overhead line of Fig. 7. Next, the cascade connection of the underground cable and overhead line is analyzed and the corresponding transient responses for different sampling rates are presented, highlighting the accuracy and effectiveness of the proposed formulation.

Fig. 8. (a) Real and (b) imaginary part of the (1,6) element of the nodal admittance matrix for the underground cable. Evaluation for different sampling rates and comparison with the corresponding direct calculation in standard NLT.

In Fig. 8 the real and imaginary part of the scalar element (1,6) of the 12 × 12 nodal admittance matrix is illustrated for the cable configuration. Results are first calculated in the frequency range of 1–100 kHz by adopting the proposed numerical improvements

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64

61

Fig. 10. Cascade connection of underground cable and overhead line.

Fig. 9. Lowest calculation accuracy of both real and imaginary part of the nodal admittance matrix for different sampling rates.

for two sampling rates, namely 3 and 10 points per decade, and are further compared to the corresponding results of Step 3 in the standard NLT where the direct calculation of all samples with frequency interval 1 Hz is performed. It is evident that the sampling rate of 10 points per decade leads to results with high accuracy, validating the proposed calculation procedure. On the other hand, the low sampling rate results in divergences from the standard NLT, which is mainly observed in the spectral peak regions. This is caused by the coarse resolution of the data at the set  in the examined frequency range. The accuracy of the proposed numerical improvements is thoroughly investigated, by evaluating the R-squared (R2 ) of (16) for different sampling rates considering the standard NLT as reference.

⎛

P/2 

⎞

⎜ (z(k) − zˆ (k)) ⎟ ⎜ ⎟ ⎜ ⎟ k=0 R2 (z) = ⎜1 − ⎟ · 100 P/2 ⎜ ⎟  ⎝ 2⎠ (z(k) − z¯ (k)) 2

The cascade connection of the short underground cable and the long overhead line implies the existence of different resonance frequencies in the transient responses. The cable coaxial propagation modes result in a resonance frequency of 44.6 kHz, indicating that a high upper frequency limit should be considered in the simulation. On the other hand, a long total observation time should be also selected, due to the line aerial propagation modes which lead to a resonance frequency of 1.5 kHz and thus to transient responses with long duration. Transient responses at ends CR1 (or Sa ) and Ra are illustrated in Figs. 11 and 12, respectively, assuming an upper frequency limit of 500 kHz and a total observation time of 15 ms. Results obtained by using the proposed numerical improvements with the sampling rate of 10 points per decade are validated with the corresponding derived by the standard NLT, where all frequency samples are directly calculated. On the other hand, the sampling rate of 3 points per decade leads to simulation errors, which result from the corresponding errors in the calculation of the cable and line nodal admittance matrices according to Fig. 9. 4.2. Computational efficiency

(16)

k=0

Here zˆ is the real or imaginary part of each interpolated element f of Y nd , z is the corresponding of the standard NLT, and z¯ is the mean of z. In Fig. 9 the lowest of the R2 measures is illustrated for both cable and overhead line configurations, as they are separately calculated for the real and imaginary part of all elements of the corresponding nodal admittance matrix. It is shown that the calculation of the nodal admittance matrix in cable and line configuration requires a sampling rate of 10 and 5 points per decade, respectively, in order to acquire a lowest R2 measure of 95%. The higher sampling rate needed in the cable configuration is due to the higher frequency dependency and scattering of the modal propagation velocities, compared to the corresponding of the overhead line. The impact of the calculation accuracy using the proposed numerical improvements is also demonstrated in a transient simulation, considering the cascade connection of the underground cable and the overhead line. As shown in Fig. 10, the cable is singlebonded at the sheath sending end and the line is left open-ended at the receiving end. For the calculation of the cable and line PUL parameters, the homogeneous earth formulations of [26,29] are used, based on simple algebraic formulas for the influence of the imperfect earth. The step voltage source applied to CS1 is considered ideal, and thus the attenuation and distortion of the responses is due to the cable and line modeling only.

In Table 1 the indicative execution times of Steps 2 and 3 are summarized for the calculation of the transient responses in Figs. 11 and 12. The remaining Steps, i.e. 1, 4 and 5, are not presented, since they are identical for both methods. All calculations are performed in MATLAB environment with a Core i7 2.93 GHz processor. The implementation of the standard NLT requires the direct calculation of the PUL parameters and propagation characteristics in 7.5·103 discrete frequencies (without counting the conjugate symmetric part) for both the underground cable and overhead line sections. This leads to a significant increase in the computational time of Step 2, being the most time-consuming step. On the other hand, the calculation of the PUL parameters and propagation characteristics in Modified Steps 2.2 and 2.3 of the proposed approach is conducted for only 51 frequencies of set  in the examined frequency range. This results in a 72.6% decrease of the total execution time. This decrease can be even higher if we consider that the number of calculated samples in the standard NLT is proportional to the total observation time and the examined frequency range, while in the proposed approach is fixed and predefined by the user in Modified Step 2.1. Next, earth formulations that involve numerical evaluation of complex integrals in TL parameter calculations are considered [30–34]. Such earth-return formulas are generally more accurate compared to simplified formulations. The relative differences on the earth-correction terms, calculated by the two formulations, are less than 1.5% below 100 kHz, but can reach up to 10% at 1 MHz [35]. Such differences may lead to divergences in the corresponding

62

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64

1.8 Standard NLT Proposed with 3 p/dec Proposed with 10 p/dec

1.6

1.4

Voltage (pu)

1.2

1

0.8

0.6

0.4

0.2

0 0

0.25

0.5

0.75

1 Time (ms)

1.25

1.5

1.75

2

Fig. 11. Transient response at end CR1 ≡ Sa , obtained from the standard NLT model and the proposed numerical improvements assuming two sampling rates. Table 1 Compared execution times considering simplified earth formulations. Standard NLT

Proposed numerical improvements (10 points per decade)

Step

Time (s)

Modified step

Time (s)

2

Cable: 21 Line: 4.5

2.2

Cable: 0.14 Line: 0.03 Cable: 4.08 Line: 1.17 Cable: 0.24 Line: 0.06 Cable: 0.98 Line: 0.77

2.3 3

3.1

Cable: 0.98 Line: 0.77

3.2 and 3.3 Total

27.25

high-frequency transient responses, especially when the induced voltages are of interest [36]. The implementation of complex integral expressions is also necessary in cases of stratified earth and in more sophisticated formulations that include admittance earth correction terms [34], since in these occasions no simple algebraic expressions exist [37]. These earth models are of major importance in electromagnetic interference (EMI) problems, where the influence of the ground and intersheath modes cannot be neglected.

Total

7.47

As an example, the step voltage transients of Fig. 10 are recalculated using the complex integral forms of the homogeneous [30] and stratified [31] earth models for the cable section, and of [32,33] for the line section, respectively. A depth of 2 m is assumed for the first layer of the stratified case, whereas an earth resistivity of 100 and 1000  m is assumed for the first and second layer, respectively. The numerical integration scheme of [21] is used for the evaluation of the PUL parameters in all earth models. The transient

Table 2 Compared execution times of the cable section considering complete earth formulations. Earth case

Homogeneous [30]

Stratified [31]

Standard NLT

Proposed numerical improvements

Step

Time (s)

Modified step

Time (s)

2

7827.94

3

0.98

Total

7828.92

2.2 2.3 3.1 3.2 and 3.3 Total

53.23 4.09 0.24 0.97 58.53

2

8733.82

3

0.98

Total

8734.8

2.2 2.3 3.1 3.2 and 3.3 Total

59.39 4.08 0.25 0.98 64.7

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64

63

The proposed methodology has a generic form so that it can be integrated with any discrete Fourier transform algorithm or any sampling scheme. Furthermore, it can be used either individually for the detailed simulation and repetitive study of TL electromagnetic transients or as a complementary tool in the validation and verification of time-domain transient models. Acknowledgments The work of A.I. Chrysochos is supported by the Alexander S. Onassis Public Benefit Foundation via a merit scholarship between 2012 and 2015. The work of T.A. Papadopoulos is supported by the ‘IKY Fellowships of Excellence for Postgraduate Studies in Greece – Siemens Program’. References Fig. 12. Transient response at end Ra , obtained from the standard NLT model and the proposed numerical improvements assuming two sampling rates.

responses are calculated with the proposed numerical improvements and the standard NLT assuming the same conditions as previously, i.e. an upper frequency limit of 500 kHz, a total observation time of 15 ms, and for the proposed method a sampling rate of 10 points per decade. In Table 2 the indicative execution times of Steps 2, 3 and Modified Steps 2.2–3.3 for the cable section only are presented for the above cases. The computational efficiency of the proposed procedure is remarkable, while the integration of these complex earth formulations with the standard NLT is practically impossible due to the unacceptable calculation burden of the PUL parameter evaluation in Step 2. 5. Conclusions In this paper, numerical improvements are proposed for the accurate and efficient frequency-domain calculation of electromagnetic transients in power TLs. The proposed procedure offers the following advantages: • The formulation of the nodal admittance matrix is performed using a minimal number of user-defined initial data, implementing a robust eigenvector tracking algorithm and a highly accurate interpolation routine. As a result, the number of required samples is not proportional to the total observation time and the examined frequency range, leading to significant reductions of the computational time in the calculation of the PUL parameters and propagation characteristics. • A minimum sampling rate of 5 and 10 points per logarithmic decade is sufficient for the accurate formulation of the nodal admittance matrix in overhead and underground configurations. In the latter case, finer sampling rates are generally required due to the stronger frequency dependency of the corresponding modal propagation characteristics. • The performance of the proposed methodology is demonstrated on a complex power line configuration, requiring both long observation time and high upper frequency limit. Results are validated with the corresponding obtained from a standard implementation of NLT, whereas the corresponding execution times are also presented. The computational efficiency of the proposed methodology is highlighted, especially in cases where complex integral formulations are used for the more accurate representation of the imperfect earth.

[1] A. Semlyen, A. Dabuleanu, Fast and accurate switching transient calculations on transmission lines with ground return using recursive convolutions, IEEE Trans. Power Appar. Syst. PAS-94 (2) (1975) 561–571. [2] J.R. Marti, Accurate modeling of frequency-dependent transmission lines in electromagnetic transient simulations, IEEE Trans. Power Appar. Syst. PAS-101 (1) (1982) 147–157. [3] L. Marti, Simulation of transients in underground cables with frequencydependent modal transformation matrices, IEEE Trans. Power Deliv. 3 (3) (1988) 1099–1110. [4] A. Morched, B. Gustavsen, M. Tartibi, A universal model for accurate calculation of electromagnetic transients on overhead lines and underground cables, IEEE Trans. Power Deliv. 14 (3) (1999) 1032–1038. [5] A. Ametani, The application of the Fast Fourier Transform to electrical transient phenomena, Int. J. Electr. Eng. Educ. 10 (4) (1973) 277–286. [6] P. Moreno, A. Ramirez, Implementation of the numerical Laplace transform: a review, IEEE Trans. Power Deliv. 23 (4) (2008) 2599–2609. [7] D.J. Wilcox, Numerical Laplace transform and inversion, Int. J. Electr. Eng. Educ. 15 (1978) 247–265. [8] S. Day, N. Mullineux, J.R. Reed, Developments in obtaining transient response using Fourier transforms. Part I: Gibbs phenomena and Fourier integrals, Int. J. Electr. Eng. Educ. 3 (1965) 501–506. [9] S. Day, N. Mullineux, J.R. Reed, Developments in obtaining transient response using Fourier transforms. Part II: use of the modified Fourier transform, Int. J. Electr. Eng. Educ. 4 (1966) 31–40. [10] P. Gómez, P. Moreno, J.L. Naredo, Frequency-domain transient analysis of nonuniform lines with incident field excitation, IEEE Trans. Power Deliv. 20 (3) (2005) 2273–2280. [11] A. Semlyen, A. Ramirez, Direct frequency domain computation of transmission line transients due to switching operations, IEEE Trans. Power Deliv. 23 (4) (2008) 2255–2261. [12] P. Gómez, F.A. Uribe, The numerical Laplace transform: an accurate technique for analyzing electromagnetic transients on power system devices, Int. J. Electr. Power Energy Syst. 31 (2009) 116–123. [13] B. Gustavsen, Validation of frequency-dependent transmission line models, IEEE Trans. Power Deliv. 20 (2) (2005) 925–933. [14] N. Nagaoka, A. Ametani, A development of a generalized frequency-domain transient program – FTP, IEEE Trans. Power Deliv. 3 (4) (1988) 1996–2004. [15] L.M. Wedepohl, C.S. Indulkar, Switching overvoltages in short crossbonded cable systems using the Fourier transform, Proc. Inst. Electr. Eng. 122 (11) (1975) 1217–1221. [16] B. Gustavsen, A.P. Brede, J.O. Tande, Multivariate analysis of transformer resonant overvoltages in power stations, IEEE Trans. Power Deliv. 26 (4) (2011) 2563–2572. [17] L.M. Wedepohl, Application of matrix methods to the solution of travellingwave phenomena in polyphase systems, Proc. Inst. Electr. Eng. 110 (12) (1963) 2200–2212. [18] J.G.F. Francis, The QR transformation, a unitary analogue to the LR transformation, Part I, Comput. J. 4 (3) (1961) 265–271. [19] J.G.F. Francis, The QR transformation, a unitary analogue to the LR transformation, Part II, Comput. J. 4 (4) (1961) 332–345. [20] A.B. Fernandes, W.L.A. Neves, Phase-domain transmission line models considering frequency-dependent transformation matrices, IEEE Trans. Power Deliv. 19 (2) (2004) 708–717. [21] G.K. Papagiannis, D.A. Tsiamitros, D.P. Labridis, P.S. Dokopoulos, Direct numerical evaluation of the earth return path impedances of underground cables, Proc. Inst. Eng. Technol. Gener. Transm. Distrib. 152 (3) (2005) 321–328. [22] F.N. Fritsch, R.E. Carlson, Monotone piecewise cubic interpolation, SIAM J. Numer. Anal. 17 (1980) 238–246. [23] D. Kahaner, C. Moler, S. Nash, Numerical Methods and Software, Prentice-Hall, Englewood Cliffs, N.J., 1989. [24] A.I. Chrysochos, T.A. Papadopoulos, G.K. Papagiannis, Robust calculation of frequency-dependent transmission line transformation matrices using

64

[25]

[26]

[27]

[28] [29]

[30] [31]

[32] [33]

[34]

[35]

A.I. Chrysochos et al. / Electric Power Systems Research 122 (2015) 56–64 the Levenberg–Marquardt method, IEEE Trans. Power Deliv. 29 (4) (2014) 1621–1629. L.M. Wedepohl, H.V. Nguyen, G.D. Irwin, Frequency-dependent transformation matrices for untransposed transmission lines using Newton–Raphson method, IEEE Trans. Power Deliv. 11 (3) (1996) 1538–1546. L.M. Wedepohl, D.J. Wilcox, Transient analysis of underground powertransmission systems: system-model and wave-propagation characteristics, Proc. Inst. Electr. Eng. 120 (2) (1973) 253–260. L.M. Wedepohl, Power system transients: errors incurred in the numerical inversion of the Laplace transform, in: Proc. 26th Midwest Symp. Circuits Systems, 1983, pp. 174–178. G. Dahlquist, A. Björk, Numerical Methods, Prentice-Hall, Englewood Cliffs, N.J., 1974, pp. 101–103. A. Deri, G. Tevan, A. Semlyen, A. Castanheira, The complex ground return plane: a simplified model for homogeneous and multi-layer earth return, IEEE Trans. Power Appar. Syst. PAS-100 (8) (1981) 3686–3693. F. Pollaczek, Über das Feld einer unendlich langen wechselstromdurchflossenen Einfachleitung, Elektr. Nachrichtentech. 3 (4) (1926) 99–139. D.A. Tsiamitros, G.K. Papagiannis, D.P. Labridis, P.S. Dokopoulos, Earth return path impedances of underground cables for the two-layer earth case, IEEE Trans. Power Deliv. 20 (3) (2005) 2174–2181. J.R. Carson, Wave propagation in overhead wires with ground return, Bell Syst. Tech. J. 5 (1926) 539–554. N. Nakagawa, A. Ametani, K. Iwamoto, Further studies on wave propagation in overhead lines with earth return: impedance of stratified earth, Proc. Inst. Electr. Eng. 120 (12) (1973) 1521–1528. T.A. Papadopoulos, D.A. Tsiamitros, G.K. Papagiannis, Earth return admittances of underground cables in non-homogeneous earth, IET Gener. Transm. Distrib. 5 (2) (2011) 161–171. O. Saad, G. Gaba, M. Giroux, A closed-form approximation for ground return impedance of underground cables, IEEE Trans. Power Deliv. 11 (3) (1996) 1536–1545.

[36] F.A. Uribe, A. Ramírez, Alternative series-based solution to approximate Pollaczek’s integral, IEEE Trans. Power Deliv. 27 (4) (2012) 2425– 2427. [37] T.A. Papadopoulos, A.I. Chrysochos, G.K. Papagiannis, Analytical study of the frequency-dependent earth conduction effects on underground power cables, IET Gener. Transm. Distrib. 7 (3) (2013) 276–287. Andreas I. Chrysochos was born in Rhodes, Greece, on October 17, 1986. He received the Dipl. Eng. degree from the School of Electrical and Computer Engineering at the Aristotle University of Thessaloniki, Greece, in 2009. Since 2010, he has been a Ph.D. candidate at the same School. His special interests are power systems modeling, computation of electromagnetic transients and power-line communications. Theofilos A. Papadopoulos was born in Thessaloniki, Greece, on March 10, 1980. He received the Dipl. Eng. and Ph.D. degrees from the School of Electrical and Computer Engineering at the Aristotle University of Thessaloniki, Greece, in 2003 and 2009 respectively. Since 2009, he has been a Post-Doc researcher at the same School. His special interests are power systems modeling, distributed generation, power-line communications and computation of electromagnetic transients. Dr. Papadopoulos received the Basil Papadias Award for the best student paper, presented at the IEEE PowerTech ‘07 Conference, Lausanne, Switzerland. Grigoris K. Papagiannis was born in Thessaloniki, Greece, on September 23, 1956. He received the Dipl. Eng. and Ph.D. degrees from the School of Electrical and Computer Engineering at the Aristotle University of Thessaloniki, Greece, in 1979 and 1998, respectively. Currently, he is an Associate Professor at the Power Systems Laboratory of the same School. His special interests are power systems modeling, computation of electromagnetic transients, distributed generation, smart grids, renewable energy sources and power-line communications.