Mathematical Problems in Engineering
Volume 2019 (2019), Article ID 4725024, 12 pages
https://doi.org/10.1155/2019/4725024

Research Article

Time Domain Electromagnetic Transient Analysis of Aerial Nonuniform Transmission Lines Excited by an Incident Electromagnetic Field

Avisaí Sánchez-Alegría, Pablo Moreno, José R. Loo-Yau, and Susana Ortega-Cisneros

Cinvestav, Unidad Guadalajara, Av. Del Bosque 1145, El Bajío, Zapopan 45019, JAL, Mexico

Correspondence should be addressed to Pablo Moreno; pmoreno@gdl.cinvestav.mx

Received 22 March 2019; Revised 12 June 2019; Accepted 9 July 2019; Published 6 August 2019

Academic Editor: Paolo Manfredi

Copyright © 2019 Avisaí Sánchez-Alegría et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In this paper a model for Nonuniform Transmission Lines for electromagnetic transient analysis that incorporates frequency dependency of electrical parameters, variation of line electrical parameters with respect to distance, and distributed excitations due to incident electromagnetic fields is presented. The model is developed using the Method of Characteristics in the actual physical domain instead of the Modal Domain; this simplifies the mathematical development and the final equations. Moreover, the equations of the resultant model are valid either for two-conductor lines or for multiconductor lines without any change; this can be an advantage when the computer programming language considers a scalar as a 1x1 matrix. The proposed numerical model is developed under the hypothesis that the dielectric surrounding the conductors is homogenous. It is shown that in this case the characteristic curves of a Nonuniform Transmission Line become straight lines. Finally, the model is validated by comparison with results obtained using the Numerical Laplace Transform method.

1. Introduction

In most studies of electromagnetic transients, the knowledge of voltages and currents at interior points of transmission lines in general has no practical importance and is not required to determine the operation state of the electrical network. However, there are some cases where it is necessary to analyze the electromagnetic behavior of the transmission line (TL) over its complete length, such as in the case of lines with distance dependent electrical parameters (nonuniform transmission lines), lines with distributed sources, lines with nonlinear electrical parameters, and lines with time dependent electrical parameters. In particular in this work the phenomena of distance dependent electrical parameters and distributed sources due to incident electromagnetic fields are considered.

Several methods have been proposed that allow computing time domain electromagnetic transients in Nonuniform Transmission Lines (NUTLs) [1–5]. In [1, 2] modeling of single-phase transmission lines was presented whereas authors in [3–5] deal with Multiconductor Transmission Lines (MTLs). In [3, 4] a two-port model is first developed in the frequency domain and afterwards it is converted to the time domain. In [5] the Finite Difference Time Domain (FDTD) Method is used to directly solve the MTL partial difference equations for the special case of frequency independent electrical parameters.

Although it is more scarce, there are some works of Nonuniform transmission lines modeling for frequency domain electromagnetic transients computation. In [6] a solution for single-phase NUTLs was presented. Later in [7] a model for external field excited Multiconductor NUTLs was developed. In such works different methods for the numerical evaluation of the inverse Laplace Transform were used.

Regarding the analysis of the effects of incident electromagnetic fields the problem can be solved directly by dividing the line in several segments and interspersing sources between line segments. The concept is simple, but it can result in a very cumbersome process for the user. Because of this several methods that aimed to model transmission lines illuminated by incident electromagnetic fields have been developed. In [8] a time domain two-port model was presented; the model is first obtained by solving the transmission line equations in the frequency domain and then a transformation to the time domain is made. In [9] the transmission line is modeled by a cascade connection of segments where each segment is represented by a macromodel. In [10] a model based on constructing a passive reduced-order system whose input/output relation approximates the Y-parameters matrix of the transmission line was presented. The model is constructed by finding an orthonormal basis for the system moments corresponding to the boundary value problem solution for the voltages and currents. The model was made independent of the shape of the external field by defining the Boundary Value Problem solution vector in terms of the transmission line state-transition matrix (or chain matrix) and the integral along the entire length of the line of the product of this matrix with an equivalent external source term. In [11–13] Finite Difference Time Domain Methods (FDTD) were developed to solve the transmission line partial differential equations. FDTD methods are versatile and relatively simple to computationally implement; however it must be taken into account that commonly they require some measures to diminish numerical oscillations.

The Method of Characteristics allows transforming the time domain transmission line partial differential equations (PDEs) to ordinary differential equations (ODEs). It has been used in the analysis of single-phase transmission lines with nonlinear capacitance due to corona effect [14], uniform MTLs with external electromagnetic field excitation [15], and nonuniform MTLs [16]. The method is simple to apply when dealing with two-conductor uniform lines, but it is a bit more complicated if multiconductor lines are analyzed because in this case modal analysis is required. The complexity increases if nonuniformities are considered because modal wave velocities become distance dependent; in a consequence the characteristics are curve lines. Moreover, modal transformation matrices are also distance dependent and therefore their derivatives with respect to distance are required [16].

In this work, the Method of Characteristics is applied to the PDEs of nonuniform Transmission Lines excited by external electromagnetic fields, where the equation of the longitudinal voltage drop is expressed in terms of the time image of the Penetration Impedance (or internal impedance), instead of the transient resistance [17], as it has been done in most of the works published up to now. Using the penetration impedance makes modal analysis unnecessary for MTLs; therefore the mathematical procedure is less cumbersome, and the resulting model is valid for two-conductor lines and for MTLs. For aerial transmission lines, which are the objective of this work, it can be considered that the dielectric that surrounds the conductors is homogeneous; in this case the characteristic curves become straight lines. This leads to a further reduction of the complexity of the numerical model but can impose a restriction for its application to other types of transmission lines.

2. Illuminated Nonuniform Transmission Line Equations

The electromagnetic behaviour of a Nonuniform MTL with frequency and distance dependent electrical parameters excited by an incident electromagnetic field, as shown in Figure 1, can be described in the Laplace Domain with the following equations [18]:

and

where and are the voltage and current vectors, respectively; , , and are the matrices of geometric inductance, capacitance, and conductance per unit length, respectively; is the penetration impedance matrix; is the distributed series voltage source that takes into account the coupling with the incident magnetic field while is a distributed shunt current source that takes into account the coupling with the incident electric field; , being ω the angular frequency and being the damping factor. For a transmission line with conductors, including the return path, matrices are of order and vectors of .

Figure 1: Illuminated Nonuniform Multiconductor Transmission line.

The geometric inductance, capacitance, and conductance matrices are defined as

where is Maxwell’s potential coefficients matrix; , , and are, respectively, the permeability, permittivity, and conductivity of the dielectric that surrounds the TL.

When the separation between conductors is large and therefore proximity effects can be neglected, the elements of P(x) are given by

where and are the radius and the height of the ith conductor, respectively; is the distance between the ith and kth conductors and is the distance between the ith conductor and the mirror image of the kth conductor. In power lines with several conductors per phase is substituted by the Geometric Mean Radius and all distances are taken from the centers of the bundles.

In accordance to Figure 1, is given by

where is the incident magnetic flux density in the z-direction corresponding to the ith conductor.

Similarly, is given as follows:

where is the incident electric field intensity in the y-direction corresponding to the ith conductor.

The penetration impedance takes into account the penetration of the electromagnetic field into the nonideal conductors of the transmission line and is given by

where is the conductor impedance matrix and the return path impedance matrix. For the case of aerial transmission lines with round conductors these matrices are calculated with well-known formulas; is calculated with formulas for cylindrical conductors in terms of modified Bessel functions [19] and using the complex penetration depth method [20]. Unfortunately, the transformation of to the time domain cannot be done analytically.

In order to facilitate the transformation to the time domain of the transmission line equations a rational approximation of can be performed using Vector Fitting (VF) [21, 22]; thus (4a) can be expressed as follows:

where is the order of the rational approximation, is the jth pole, and is the corresponding residues matrix; is a real constant matrix.

Passivity enforcement of models synthesized in the frequency domain with rational functions is of paramount importance when simulating electromagnetic transients in the time domain, since it has been found that when the fitted model is not passive, even when the model contains stable poles only, numerical instabilities arise [23, 24]. Passivity requires the real part of the fitted at all frequencies to be Positive Definite. In the VF method first a rational approximation with stable poles is obtained; then a correction is calculated that enforces the Positive Definite criterion to be satisfied; the enforcement procedure is based on linearization and constrained minimization by Quadratic Programming [24]. It should be mentioned that the elements of the penetration impedance matrix of aerial power transmission lines are smooth functions of the frequency and that the authors have not yet encountered problems of passivity in the rational expansions.

From (4b) it can be observed that the direct current resistance of the transmission line is given by

When performing the rational approximation, the condition (4c) is very difficult to achieve, especially when the maximum frequency is high. Enforcement of (4c) requires an extremely low value of the initial frequency for the rational approximation; this causes an increase of the error at high frequencies. An approach that represents a compromise between low and high frequency errors uses a not very low value of the initial frequency and then enforces (4c) by correcting the value of .

Substituting (4b) into (1a) yields

Transforming (5a) to the time domain the following can be written:

where

where represents the inverse Laplace Transformation

Expression (6b) represents a time domain convolution that can be computed with a recursive algorithm; thus (5b) becomes

where

The shunt conductance is negligible in aerial transmission lines; for this case transforming (1b) to the time domain gives

where

3. Method of Characteristics: A Nonmodal Approach

Expressions (7a) and (8a) can be grouped as follows:

where

Matrices A and B are of order 2N×2N and vectors W, H, and F of 2N×1.

System (9a) is hyperbolic if A is diagonalizable with real eigenvalues [25]; in other words there is a matrix function E(x) such that

where are real and the norms of and are bounded in x, for .

It can be shown that the matrix of eigenvectors is given by

where U is the N×N identity matrix and R₀(x) is the Characteristic Resistance of the line defined as follows:

In (11b) is the intrinsic impedance of the dielectric that surrounds the conductors given by

In addition, it can be proved that the eigenvalues of A are given by

where is a diagonal matrix given by

where is the velocity of the wave in the dielectric defined as follows:

In accordance with (12b) and (12c) for the positive eigenvalues in (12a) the following equation can be written:

and for the negative eigenvalues of (12a)

where

Solutions to (13a) and (13b) are the Characteristics of the hyperbolic system (9a). Since is diagonal with all its elements equal, (13a) and (13b) are actually scalar equations that can be expressed as follows:

and

Solving (14a) and (14b) gives

and

where α and β are real integration constants and

Equations (14c) and (14d) represent two families of Characteristic curves in the x-t plane, as shown in Figure 2. These characteristics have been obtained from the product , and since it is diagonal it was not necessary to obtain its eigenvalues; consequently it is not necessary to define modal voltages and currents.

Figure 2: Characteristic curves in the x-t plane for a NUTL.

On a previous work of the application of the method of characteristics to solve the equations of Nonuniform MTLs [16], the starting point was the Transmission Line Equations in terms of the Transient Resistance; such approach requires using Modal Analysis, which yields N different Characteristics with positive slope and N different Characteristics with negative slope.

In the case of a nonuniform transmission line where the dielectric that surrounds the conductors is homogeneous, i.e., ε and μ are constants, it is clear from (12c) that the eigenvalues do not depend on x and the solutions of (13a) and (13b) reduce to the solutions of a Uniform Transmission line as follows:

and

Therefore, for this special case the characteristics are straight lines regardless of the variations of heights or radii of conductors or distance between them.

Left multiplying (9a) by E(x) gives

and

In (16a) and (16b) it has been taken into account that

and

Considering definitions (13a) and (13b) the terms in parenthesis in (16a) and (16b) become total derivatives; thus the following can be written:

and

Expressions (19a) and (19b) along with (13a) and (13b) are ODEs that are equivalent to the PDEs system (9a). In (19a) and (19b) it can be observed that the Nonmodal approaches and , which take into account the effects of the penetration of the electromagnetic field in nonideal conductors, are out of the differential operations; therefore they do not change the main features of the propagation; in other words, their effects can be seen as a distortion of the wave that travels along a lossless nonuniform transmission line.

Different from the method presented in [16] (19a) and (19b) are in terms of actual voltages and currents; moreover there are no modal transformation matrices inside differential operations, as in of [16]. This greatly simplifies the numerical solution.

4. Numerical Solution

4.1. Interior Points

Numerical solution of the line model given by (19a) and (19b) is performed using finite differences. In the general case of a Nonuniform Line the Characteristics are curve lines and strictly a nonuniform mesh is required. A more practical option is to define a uniform mesh and then perform interpolations to determine the quantities at the required crossing points for the numerical method; this is shown in Figure 3 where the values of voltages and currents at points B’ and C’ should be calculated with the known values at points B, D, and C. As discussed previously, in the special case of an aerial transmission line the Characteristics are straight lines and then interpolations are not required.

Figure 3: Discretization mesh.

Consider that solutions to (19a) and (19b) are known at time instant (points B, C, and D) and that these solutions are to be extended to time instant t. Note that the relation between ∆t and ∆x must satisfy the Courant-Friedrichs-Lewy (CFL) condition [25]. Using a central differences scheme on (19a) and (19b), the following is obtained:

and

where

where subscripts A, B, and C indicate that calculations are performed using the transmission line parameters at that point. For the incident electromagnetic field sources, subscripts AB and AC indicate that calculation is made at the center of the corresponding characteristic.

Matrix equations (20a) and (20b) can be grouped as follows:

where

By solving (22a) and can be found. Sources and take into account the effect of the incident electromagnetic field, and they are calculated in accordance with the orientation of the components of the incident electromagnetic field or when this field is not present they are simply given a value of zero.

4.2. End Points

At each end of the transmission line there exists only one characteristic, as shown in Figure 4; therefore one more expression for each line end is needed. In the case of a radial transmission system simple circuits can be connected at the source and the load, as shown in Figure 5. In accordance with Figure 5 for the sending end the following can be written:

and for the receiving end

(a)

(b)

Figure 4: Conditions at Transmission Line Ends: (a) sending end and (b) receiving end.

(a)

(b)

Figure 5: Simple terminations of a transmission line: (a) sending end and (b) receiving end.

To obtain the voltage and current at the sending end (23a) is solved simultaneously with the second row of (22a), as shown in (24a).

and for the receiving end (23b) is solved simultaneously with the first row of (22a) as follows:

When the objective is the study of a transmission network the method described above is not practical. A better approach is to develop a two-port model for the transmission line ends that can be included in programs like ATP/EMTP or SPICE. From the second row of (22a) the following can be written for the sending end:

where

Similarly, from the first row of (22a) the following can be written for the receiving end:

where

Expressions (25a) and (26a) represent Norton models for the transmission line ends, as presented in Figure 6. In (25a) and (26a) the currents and are defined as positive when they enter the transmission line, as is also shown in Figure 6.

Figure 6: Norton equivalent circuit for a transmission line.

Using the model presented in Figure 6 into a network requires two stages: given initial or previous network conditions, solve for the currents and voltages of the network by means of the Nodal Method or the Modified Nodal Method and calculate currents and voltages at interior point of the transmission line (or lines); then the network can be solved again for the next time step. Regarding the second stage of the procedure, only the first interior point after the sending end and the last interior point before the receiving end is actually required to proceed to the following network solution stage. Therefore, to speed the calculation and when the adequate tools are available the rest of the interior points can be computed in parallel to the network solution.

5. Study Cases

To verify the model developed in this work, a comparison with results obtained using the Numerical Laplace Transform [7, 26] was performed. Three study cases are presented; the first one is a single-phase transmission line with catenary; the second one is a MTL crossing over a river and the last one is a MTL illuminated by an incident electromagnetic field. In the three study cases the Penetration Impedance was synthetized with VF using 100 points logarithmically spaced in the frequency band 1 Hz - 6 MHz

5.1. Nonuniform Single-Phase Transmission Line

A single-phase line long with sagging between towers as shown in Figure 7 is considered. The conductor is made of aluminum with a radius of , the line maximum height is at the towers, and the minimum height is at the middle of the span; a ground resistivity of is assumed. A voltage with waveform of double linear ramp is applied at the sending node. The source impedance is and a load of is connected at the line receiving end. For the simulation the number of samples was 1024 and the time step was . In this case the synthesis with VF yielded 10 poles.

Figure 7: Sagging between towers.

Figure 8 shows the voltages at the sending and receiving nodes; it is observed that the waveforms obtained with the Method of Characteristics (MC) and the Numerical Laplace Transform (NLT) are practically the same.

Figure 8: Voltages at the sending and receiving ends of the monophasic nonuniform line.

5.2. River Crossing Nonuniform MTL

As a second example, consider a nonuniform horizontal three-phase TL crossing a river, as shown in Figure 9 [7]. The line is long and has one conductor per phase; the radius of the conductors is and the phase separation is ; a water resistivity of is assumed. The sending end is at a height of and the receiving end is at . The voltage source is a unit step applied to all conductors at the sending node. A very small source impedance is considered () and the receiving end is left open.

Figure 9: Transmission line crossing a river.

The number of samples was 1024 with a time step of . To fit the Penetration Impedance VF delivered 18 poles. Figures 10 and 11 present the waveforms of phases A (external conductor) and B (central conductor), respectively, at the receiving end obtained with the MC and the NLT; it is observed again that both methods produce very close results.

Figure 10: Voltages at phase A of the receiving end for the transmission line crossing a river.

Figure 11: Voltages at phase B of the receiving end for the transmission line crossing a river.

5.3. Nonuniform MTL Excited by an Incident Electromagnetic Field

Consider a horizontal three-phase transmission line segment with catenary, as shown in Figure 12. The line is long and has one conductor per phase; the phase separation is and the radius of the conductors is ; it is assumed that the earth resistivity is . A load of is connected to each conductor at the transmission line ends. For this example, the penetration impedance was synthesized with 18 poles.

Figure 12: Transmission line segment excited by an incident electromagnetic field.

For this study case a plane wave propagating in the x-direction was considered, composed by an electric field oriented in the y-direction given by

and a magnetic field oriented in the z-direction was expressed as follows:

where and .

Figures 13–15 show induced voltages at the center conductor. In Figure 13 the voltage induced on the load resistances at the sending and the receiving ends is shown due to the incident electromagnetic field. Figure 14 shows the induced voltages due to only the electric field and Figure 15 the voltages induced by the magnetic field.

Figure 13: Induced voltages at the transmission line ends due to the incident electromagnetic field.

Figure 14: Induced voltages at the transmission line ends due to the incident electric field.

Figure 15: Induced voltages at the transmission line ends due to the incident magnetic field.

6. Conclusions

In this work a model for aerial NUTLs illuminated by incident electromagnetic fields for time domain electromagnetic transient studies has been presented. The model was developed using an alternative approach of the Method of Characteristics; this approach takes as starting point the longitudinal telegrapher’s equation expressed in terms of the penetration impedance, instead of the transient resistance widely used in previous works.

With the approach of this work the characteristics of the transmission line PDEs are obtained from the product , which is diagonal; thus there is no need of using modal analysis. It has been shown that when the dielectric that surrounds the TL conductors is homogeneous the characteristics are straight lines as in the case of uniform TLs; the characteristics become curve lines when the permittivity and the permeability of the dielectric are distance dependent or at least one of them.

From the development presented it can be said that the presence of nonideal conductors does not change the main features of the propagation phenomenon and their effects appear as a distortion of the traveling wave. It can also be said that when the method of characteristics based on modal analysis is used, the appearance of several modes of propagation is due to the need to take into account the distortion of the traveling wave produced by the penetration of the electromagnetic field in nonideal conductors.

The proposed model retains the low complexity typical of FDTD methods with the advantage that until now it has been shown that the Method of Characteristics is less prone to numerical oscillations. Moreover, the model allows knowing the voltage and current profiles along the whole transmission line, which can be important in insulation coordination studies of power transmission lines and probably of other power apparatuses modeled with distributed parameters.

Data Availability

Data used to support the findings of this work are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

H. Nguyen, H. Dommel, and J. Marti, “Modelling of single-phase nonuniform transmission lines in electromagnetic transient simulations,” IEEE Transactions on Power Delivery, vol. 12, no. 2, pp. 916–921, 1997.

M. de Barros and M. Almeida, “Computation of electromagnetic transients on nonuniform transmission lines,” IEEE Transactions on Power Delivery, vol. 11, no. 2, pp. 1082–1091, 1996.

A. Ramirez, A. Semlyen, and R. Iravani, “Modeling nonuniform transmission lines for time domain simulation of electromagnetic transients,” IEEE Transactions on Power Delivery, vol. 18, no. 3, pp. 968–974, 2003.

A. C. S. Lima, R. A. R. Moura, B. Gustavsen, and M. A. O. Schroeder, “Modelling of non-uniform lines using rational approximation and mode revealing transformation,” IET Generation, Transmission & Distribution, vol. 11, no. 8, pp. 2050–2055, 2017.

K. Afrooz and A. Abdipour, “Efficient method for time-domain analysis of lossy nonuniform multiconductor transmission line driven by a modulated signal using FDTD technique,” IEEE Transactions on Electromagnetic Compatibility, vol. 54, no. 2, pp. 482–494, 2012.

E. Oufi, A. Alfuhaid, and M. Saied, “Transient analysis of lossless single-phase nonuniform transmission lines,” IEEE Transactions on Power Delivery, vol. 9, no. 3, pp. 1694–1700, 1994.

P. Gomez, P. Moreno, and J. Naredo, “Frequency-domain transient analysis of nonuniform lines with incident field excitation,” IEEE Transactions on Power Delivery, vol. 20, no. 3, pp. 2273–2280, 2005.

A. Andreotti, A. Pierno, and V. A. Rakov, “A new tool for calculation of lightning-induced voltages in power systems—part i: development of circuit model,” IEEE Transactions on Power Delivery, vol. 30, no. 1, pp. 326–333, 2015.

X. Liu, X. Cui, and L. Qi, “Calculation of lightning-induced overvoltages on overhead lines based on DEPACT macromodel using circuit simulation software,” IEEE Transactions on Electromagnetic Compatibility, vol. 54, no. 4, pp. 837–849, 2012.

E. Gad, “Circuit-based analysis of electromagnetic field coupling with nonuniform transmission lines,” IEEE Transactions on Electromagnetic Compatibility, vol. 50, no. 1, pp. 149–165, 2008.

M. Brignone, F. Delfino, R. Procopio, M. Rossi, and F. Rachidi, “Evaluation of Power system lightning performance, part i: model and numerical solution using the PSCAD-EMTDC platform,” IEEE Transactions on Electromagnetic Compatibility, vol. 59, no. 1, pp. 137–145, 2017.

S. I. Abouzeid, G. Shabib, and A. Z. E. D. Mohamed, “Induced voltages on overhead transmission lines because of nearby included lightning channel,” IET Generation, Transmission & Distribution, vol. 9, no. 13, pp. 1672–1680, 2015.

E. L. Tan and Z. Yang, “Non-uniform time-step FLOD-FDTD method for multiconductor transmission lines including lumped elements,” IEEE Transactions on Electromagnetic Compatibility, vol. 59, no. 6, pp. 1983–1992, 2017.

A. Ramírez, J. L. Naredo, P. Moreno, and L. Guardado, “Electromagnetic transients in overhead lines considering frequency dependence and corona effect via the method of characteristics,” International Journal of Electrical Power & Energy Systems, vol. 23, no. 3, pp. 179–188, 2001.

A. Ramirez, J. Naredo, and P. Moreno, “Full frequency-dependent line model for electromagnetic transient simulation including lumped and distributed sources,” IEEE Transactions on Power Delivery, vol. 20, no. 1, pp. 292–299, 2005.

A. Chavez, P. Moreno, and J. Naredo, “Fast transient analysis of nonuniform multiconductor frequency-dependent transmission lines,” IEEE Transactions on Power Delivery, vol. 7, no. 11, pp. 1185–1193, 2013.

R. Radulet, A. L. Timotin, and A. Tugulea, “The propagation equations with transient parameters for long lines with losses,” Revue Roumaine des Sciences Techniques-Serie Electrotechnique et Energetique, vol. 15, no. 4, pp. 587–599, 1970.

C. D. Taylor, R. S. Satterwhite, and C. W. Harrison, “The response of a terminated two-wire transmission line excited by a nonuniform electromagnetic field,” IEEE Transactions on Antennas and Propagation, vol. 13, no. 6, pp. 987–989, 1965.

N. D. Tleis, Power Systems Modelling and Fault Analysis, Theory and Practice, Newnes-Elsevier, Oxford, UK, 2008.

J. A. Martinez-Velasco, A. I. Ramirez, and M. Dávila, Overhead Lines, Power System Transients: Parameter Determination, CRC Press-Taylor & Francis Group, Boca Raton, Fla, USA, 2010.

B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Transactions on Power Delivery, vol. 14, no. 3, pp. 1052–1061, 1999.

“The Vector Fitting,” https://www.sintef.no/projectweb/vectfit/.

A. Dounavis, R. Achar, and M. S. Nakhla, “Efficient passive circuit models for distributed networks with frequency-dependent parameters,” IEEE Transactions on Advanced Packaging, vol. 23, no. 3, pp. 382–392, 2000.

B. Gustavsen and A. Semlyen, “Enforcing passivity for admittance matrices approximated by rational functions,” IEEE Transactions on Power Systems, vol. 16, no. 1, pp. 97–104, 2001.

J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Wadsworth & Brooks, Pacific Grove, Calif, USA, 1989.

P. Moreno and A. Ramírez, “Implementation of the numerical laplace transform: a review,” IEE Trans. Power Delivery, vol. 23, no. 4, pp. 2599–2609, 2008.