Computation of Transient Profiles along Nonuniform Transmission Lines Including Time-Varying and Nonlinear Elements Using the Numerical Laplace Transform †

Electromagnetic transients are responsible for overvoltages and overcurrents that can have a negative impact on the insulating elements of the electrical transmission system. In order to reduce the damage caused by these phenomena, it is essential to accurately simulate the effect of transients along transmission lines. Nonuniformities of transmission line parameters can affect the magnitude of voltage transients, thus it is important to include such nonuniformities correctly. In this paper, a frequency domain method to compute transient voltage and current profiles along nonuniform multiconductor transmission lines is described, including the effect of time-varying and nonlinear elements. The model described here utilizes the cascade connection of chain matrices in order to take into consideration the nonuniformities along the line. This technique incorporates the change of parameters along the line by subdividing the transmission line into several line segments, where each one can have different electrical parameters. The proposed method can include the effect of time-dependent elements by means of the principle of superposition. The numerical Laplace transform is applied to the frequency-domain solution in order to transform it to the corresponding time-domain response. The results obtained with the proposed method were validated by means of comparisons with results computed with ATP (Alternative Transients Program) simulations, presenting a high level of agreement.


Introduction
Electromagnetic transients can produce overvoltages and overcurrents that can have a negative impact on electric power systems. In order to reduce the potential deterioration or damage due to this condition, accurate transient simulations are needed [1]. Typically, the transient analysis of electrical systems is performed by means of two-port models of uniform transmission lines, thus the voltage/current measurements are available at certain nodes of the network. However, the maximum transient overvoltages and overcurrents may appear at interior points of the transmission line [2]. In such cases, the traditional simulation methods may not be well suited to correctly analyze this kind of phenomena. Additionally, nonuniformities can be present along the transmission line in the form of parameter variations such as the height of the line or the properties of the terrain; if these nonuniformities are too prominent, the distribution and magnitude of the voltage and currents can drastically be affected along the transmission line in comparison with uniform transmission lines (lines with space-independent per-unit-length parameters) [3,4].
There has been a considerable interest in developing methods to accurately model nonuniform transmission lines during electromagnetic transients. Previous works have presented line models applying different techniques such as the numerical Laplace transform, the method of characteristics, rational approximations, among others, with good results [5][6][7][8][9][10][11][12]. However, in general, these methods are only able to provide voltage and current information at the ends of the line, which in some cases may not be enough to correctly analyze the transient behavior of a transmission line [2], such as for insulation design or protection purposes.
A frequency domain method for the computation of transient voltage and current profiles along transmission lines was reported in [13,14], and later extended to include time-varying and nonlinear elements [15]. However, this method cannot be applied to nonuniform transmission lines. In [16], Laplace-domain and time-domain methods were tested in the computation of transient profiles on a nonuniform electronic system. It was found that the line model that used the inverse numerical Laplace transform (INLT) provided the most accurate results. However, the method described in [16] can only be applied to time-invariant linear systems.
Expanding upon the aforementioned publications, the main contribution of the present paper is the complete description and verification of a method for the computation of transient voltage and current profiles along nonuniform transmission lines. This method utilizes a modeling approach defined in the frequency domain and based on the cascaded connection of chain matrices. Furthermore, using the superposition technique, the proposed method can include time-dependent and nonlinear elements (such as switching devices and surge arresters) in the computation of transient profiles along nonuniform lines, something that has not been conducted in any previously reported research. Since the line model is defined in the frequency domain, it can take into account the frequency dependence of electrical parameters in a straightforward manner, providing more accurate results in comparison with existing methods defined in the time domain.
The method presented here makes use of the inverse numerical Laplace transform [17,18] to convert the computed frequency domain solution to a time-domain transient response. This method has strong potential for application in fault location and insulation coordination, with particular accuracy benefits for lines with prominent nonuniformities, such as river crossings, hilly terrains, and other substantial sagging conditions, which are commonly encountered in large countries such as China, Canada, India, Russia, and Brazil. For example, very challenging river crossings are found in Brazil for overhead lines constructed to connect the power generation in the Amazon Basin to the main load centers of the country. These river crossings are in the order of 2 km leading to very tall towers and wide line spans [19]. Accurate fault location under these circumstances requires an appropriate consideration of wave propagation along nonuniform lines, which can be achieved with the method described here. Additionally, the proposed method can be expanded to the modeling of other power system nonuniform elements, such as transmission towers [20] and rotating machines [21].
In order to validate the accuracy of the proposed method, the results obtained are compared with those obtained from simulations performed with ATP (Alternative Transients Program) [22]. In the ATP simulations, the J. Marti line model [23] was used, and the transmission line was subdivided into several line segments to allow the connection of measuring probes at internal points of the line, as well the inclusion of nonuniformities.
The computation of transient profiles along nonuniform transmission lines including nonlinear and time-varying conditions has not been previously reported, providing an original contribution to the current state of the art of the topic.

Nonuniform Transmission Line Model for the Transients Profiles Computation
This section describes the transmission line model used to compute the transient profiles as well as the technique used to include time-varying and nonlinear elements. Additionally, a brief explanation for the implementation of the INLT algorithm is presented.

Nonuniform Transmission Line Model
This work introduces the nonuniformities along transmission lines by means of the technique of cascade connection of chain matrices, as it has been previously shown to be an effective technique for the simulation of electromagnetic transients when nonuniform transmission lines are considered [11,24].
Initially, the uniform transmission line of Figure 1 is considered. This line can be represented by a two-port model known as transfer or ABCD matrix model: where In Equation (2) Z 0 , Y 0 , Z, Y, and L are the line's characteristic impedance, characteristic admittance, series impedance, shunt admittance and length, respectively. Additionally, the Laplace variable is defined as s = c + jω, where ω is given by 2π f .

Nonuniform Transmission Line Model for the Transients Profiles Computation
This section describes the transmission line model used to compute the transient profiles as well as the technique used to include time-varying and nonlinear elements. Additionally, a brief explanation for the implementation of the INLT algorithm is presented.

Nonuniform Transmission Line Model
This work introduces the nonuniformities along transmission lines by means of the technique of cascade connection of chain matrices, as it has been previously shown to be an effective technique for the simulation of electromagnetic transients when nonuniform transmission lines are considered [11,24].
Initially, the uniform transmission line of Figure 1 is considered. This line can be represented by a two-port model known as transfer or ABCD matrix model: In Equation (2) Z0, Y0, Z, Y, and L are the line's characteristic impedance, characteristic admittance, series impedance, shunt admittance and length, respectively. Additionally, the Laplace variable is defined as = + , where is given by 2 . By changing the direction of IL in Figure 1, Equation (1) is modified in the following manner: or in a compact form: Matrix Φ in Equation (4) is the chain matrix of the transmission line. Due to the fact that the currents at both ends of the line have the same direction, multiple transmission lines can be cascade-connected by using their corresponding chain matrices, as it can be seen in Figure 2. By changing the direction of I L in Figure 1, Equation (1) is modified in the following manner: or in a compact form: Matrix Φ in Equation (4) is the chain matrix of the transmission line. Due to the fact that the currents at both ends of the line have the same direction, multiple transmission lines can be cascade-connected by using their corresponding chain matrices, as it can be seen in Figure 2.  Figure 2 can also be interpreted as one transmission line subdivided into several smaller line segments, where each segment is represented by a unique chain matrix with its own independent electrical properties. Using this approach, it is possible to build a transmission line model that includes nonuniformities. With this consideration in mind, it can be deduced from Figure 2 that the voltage and current at the beginning of each line segment and the chain matrix of the same segment can be used to compute the voltage and current of the next segment: or in a general way: Equation (6) can be used to compute the voltage and current profiles along a nonuniform transmission line. The transient profiles are computed in a sequential manner, obtaining the voltage and current at the end of the first line segment from the chain matrix and from the voltages and currents at the beginning of the same segment; this process is repeated until the voltage and current along the whole line have been computed. It can also be observed in Equation (6) that V0 and I0 are required as initial values of the algorithm; in order to compute such values, a two-port admittance representation of the complete transmission line from the chain matrix ΦFL is defined as:    Figure 2 can also be interpreted as one transmission line subdivided into several smaller line segments, where each segment is represented by a unique chain matrix with its own independent electrical properties. Using this approach, it is possible to build a transmission line model that includes nonuniformities. With this consideration in mind, it can be deduced from Figure 2 that the voltage and current at the beginning of each line segment and the chain matrix of the same segment can be used to compute the voltage and current of the next segment: or in a general way: Equation (6) can be used to compute the voltage and current profiles along a nonuniform transmission line. The transient profiles are computed in a sequential manner, obtaining the voltage and current at the end of the first line segment from the chain matrix and from the voltages and currents at the beginning of the same segment; this process is repeated until the voltage and current along the whole line have been computed. It can also be observed in Equation (6) that V 0 and I 0 are required as initial values of the algorithm; in order to compute such values, a two-port admittance representation of the complete transmission line from the chain matrix Φ FL is defined as: where V 0 (s) is obtained by solving (8) for the voltages vector. I 0 (s) is computed as follows:

Modeling of Time-Varying Elements
It can be difficult to include time-varying conditions, such as switching maneuvers, when methods defined in the frequency domain are used to simulate electromagnetic transients. However, the principle of superposition has demonstrated to be an efficient method to include such conditions in the frequency domain [25], as explained below.
The closing of a switch can be computed using the circuits presented in Figure 3, which represent the state of the circuit before and after the closing maneuver. V E , V S , and V R are the voltage at the source side, at the line's sending node, and at the line's receiving node, respectively. The nodal voltage vector V NL (s) is formed by the three subvectors in Figure 3 as shown below: and can be obtained from the following expression [11]: Energies 2019, 12, x FOR PEER REVIEW 5 of 16

Modeling of Time-Varying Elements
It can be difficult to include time-varying conditions, such as switching maneuvers, when methods defined in the frequency domain are used to simulate electromagnetic transients. However, the principle of superposition has demonstrated to be an efficient method to include such conditions in the frequency domain [25], as explained below.
The closing of a switch can be computed using the circuits presented in Figure 3, which represent the state of the circuit before and after the closing maneuver. VE, VS, and VR are the voltage at the source side, at the line's sending node, and at the line's receiving node, respectively. The nodal voltage vector VNL(s) is formed by the three subvectors in Figure 3 as shown below: and can be obtained from the following expression [11]: In Equation (12), Ybus0 is the nodal admittance matrix before the switch closing (Figure 3a), Ybus1 is the admittance matrix modified by the closing operation (Figure 3b), IN0 contains the initially injected currents, and IN1 is the injection current vector due to the switch operation. An example of the construction of the nodal admittance matrices Ybus0 and Ybus1 for the transmission line in Figure 3 is presented below: If Ybus0 (Figure 3a) is to be built using (13), the elements , , , and would assume the following values: where is the open switch's admittance between its poles, ideally = 0. = + , with being the admittance connected at the sending node, and = + where represents the self-admittance of the transmission line. The matrix Ybus1 (Figure 3b) can be built in a similar way, but replacing by , that is, the admittance of the switch when closed.
The procedure presented above can be extended to any number of changes in the circuit topology using the following expression: where Ybusi and INi are the modified admittance matrix and the injection current vector corresponding to the n-th switch operation, respectively. A comprehensive explanation of this procedure can be found in [11]. V0(s) in (11) is the voltage at the sending node of the line needed in (6) to compute the transient voltage and current profiles. Meanwhile, the current at the beginning of the line is computed with (10) using V0(s) and VL(s) from (11). In Equation (12), Ybus 0 is the nodal admittance matrix before the switch closing (Figure 3a), Ybus 1 is the admittance matrix modified by the closing operation (Figure 3b), I N0 contains the initially injected currents, and I N1 is the injection current vector due to the switch operation. An example of the construction of the nodal admittance matrices Ybus 0 and Ybus 1 for the transmission line in Figure 3 is presented below: If Ybus 0 (Figure 3a) is to be built using (13), the elements Y 11 , Y 12 , Y 21 , and Y 22 would assume the following values: with Y s being the admittance connected at the sending node, and Y 22 = Y LL + Y switcho where Y LL represents the self-admittance of the transmission line. The matrix Ybus 1 (Figure 3b) can be built in a similar way, but replacing Y switcho by Y switchc , that is, the admittance of the switch when closed.
The procedure presented above can be extended to any number of changes in the circuit topology using the following expression: where Ybus i and I Ni are the modified admittance matrix and the injection current vector corresponding to the n-th switch operation, respectively. A comprehensive explanation of this procedure can be found in [11].
V 0 (s) in (11) is the voltage at the sending node of the line needed in (6) to compute the transient voltage and current profiles. Meanwhile, the current at the beginning of the line is computed with (10) using V 0 (s) and V L (s) from (11).

Inclusion of Nonlinear Elements
It has been demonstrated that the principle of superposition can be used to overcome the difficulties of the inclusion of nonlinear components in the frequency domain [25]. This is achieved by means of a sequence of switching operations.
In order to include a nonlinear element in a simulation performed using a method based in the frequency domain, first, its nonlinear characteristic curve must be approximated by a piece-wise curve, as illustrated in Figure 4.

Inclusion of Nonlinear Elements
It has been demonstrated that the principle of superposition can be used to overcome the difficulties of the inclusion of nonlinear components in the frequency domain [25]. This is achieved by means of a sequence of switching operations.
In order to include a nonlinear element in a simulation performed using a method based in the frequency domain, first, its nonlinear characteristic curve must be approximated by a piece-wise curve, as illustrated in Figure 4.
By approximating the v-i characteristic of a nonlinear element as shown above, such an element can be represented as a network of N parallel-connected branches, as shown in Figure 5. In the circuit presented in Figure 5, the switch in the n-th branch operates according to the reference voltage; it closes when the voltage across nodes j and k goes above Vrefn and it opens when the voltage drops below Vrefn. It is important to mention that the switches must operate in a successive manner, meaning that a switch can only close when the one in the preceding branch is closed, and a switch can only open when the one in the following branch is open.
When the switch connected to the n-th branch is closed, the equivalent Thevenin resistance of the circuit in Figure 5 must be equal to the value of , that is, the slope of the n-th linear segment in The element R n of the n-th linear element from Figure 4 and the voltage V n are computed as follows: By approximating the v-i characteristic of a nonlinear element as shown above, such an element can be represented as a network of N parallel-connected branches, as shown in Figure 5.

Inclusion of Nonlinear Elements
It has been demonstrated that the principle of superposition can be used to overcome the difficulties of the inclusion of nonlinear components in the frequency domain [25]. This is achieved by means of a sequence of switching operations.
In order to include a nonlinear element in a simulation performed using a method based in the frequency domain, first, its nonlinear characteristic curve must be approximated by a piece-wise curve, as illustrated in Figure 4.
By approximating the v-i characteristic of a nonlinear element as shown above, such an element can be represented as a network of N parallel-connected branches, as shown in Figure 5. In the circuit presented in Figure 5, the switch in the n-th branch operates according to the reference voltage; it closes when the voltage across nodes j and k goes above Vrefn and it opens when the voltage drops below Vrefn. It is important to mention that the switches must operate in a successive manner, meaning that a switch can only close when the one in the preceding branch is closed, and a switch can only open when the one in the following branch is open.
When the switch connected to the n-th branch is closed, the equivalent Thevenin resistance of the circuit in Figure 5 must be equal to the value of , that is, the slope of the n-th linear segment in In the circuit presented in Figure 5, the switch in the n-th branch operates according to the reference voltage; it closes when the voltage across nodes j and k goes above V refn and it opens when the voltage drops below V refn . It is important to mention that the switches must operate in a successive manner, meaning that a switch can only close when the one in the preceding branch is closed, and a switch can only open when the one in the following branch is open.
When the switch connected to the n-th branch is closed, the equivalent Thevenin resistance of the circuit in Figure 5 must be equal to the value of R n , that is, the slope of the n-th linear segment in the approximation of Figure 4: Finally, by solving (17) for R xn , the value of such resistance can be computed: V xn is computed as: With this approach, it is possible to compute the transient profiles along nonuniform transmission lines with nonlinear and/or time-varying elements using (6), calculating the voltage and current at the sending node with (10) and (14), and implementing the nonlinear element model presented in this section.

Inverse Numerical Laplace Transform
As the last step, the INLT is applied to the computed transient profiles in the frequency domain with (6); this is done in order to transform the results to the time domain. This method has proven to be very accurate for the study of electromagnetic transients [11,24,25]. The implementation of the INLT algorithm is described below in a brief manner; references [17,18] can be consulted for a thorough explanation. Considering a time function f (t) as a real and causal function, the inverse Laplace transform can be written as: In this work, the numerical evaluation of (20) is done considering an odd sampling in the frequency spectrum (using a spacing of 2∆ω), and normal time steps ∆t in the time domain. With these considerations in mind, the following definitions for the discrete functions in the time and frequency domain for N equally spaced samples are made: where n, m = 0, 1, 2, . . . , N − 1 and chat reduces the aliasing errors of the algorithm; in this work, c is defined as c = 2∆ω. Additionally, it is necessary to establish finite integration limits for the numerical evaluation of (20), that is, the maximum frequency Ω and the observation time T. The observation time is computed from: and the following relations can be established: Finally, by implementing the odd sampling presented in (21) and including a window function σ m , the Laplace transform in (20) can be numerically approximated by: for n, m = 0, 1, 2, . . . , N -1 (24) where: In Equation (24), the term inside the brackets corresponds to the fast Fourier transform algorithm; this allows computing time savings if N is equal to an integer power of two.

Test Cases
Three test cases are presented in order to validate the proposed method. The frequency dependence of the line's electrical parameters was taken into account by means of the application of the complex image method to introduce the earth-return impedance, as well as the complex penetration depth to include the skin effect in the line conductors [26]. The accuracy of the method was validated through comparisons with measurements from simulations performed with ATP, where the frequency-dependent J. Marti line model was discretized in several smaller segments to allow the introduction of nonuniformities and the connection of measuring probes at interior points. The length of the simulated lines in the test cases was kept short in order to reduce the implementation burden in ATP; with the proposed method, the length and discretization of the line did not represent a problem.

Highly Nonuniform Transmission Line
A highly nonuniform three-phase transmission line was considered. The nonuniformity was included by the sagging of the conductors as shown in Figure 6.
In Equation (24), the term inside the brackets corresponds to the fast Fourier transform algorithm; this allows computing time savings if is equal to an integer power of two.

Test Cases
Three test cases are presented in order to validate the proposed method. The frequency dependence of the line's electrical parameters was taken into account by means of the application of the complex image method to introduce the earth-return impedance, as well as the complex penetration depth to include the skin effect in the line conductors [26]. The accuracy of the method was validated through comparisons with measurements from simulations performed with ATP, where the frequency-dependent J. Marti line model was discretized in several smaller segments to allow the introduction of nonuniformities and the connection of measuring probes at interior points. The length of the simulated lines in the test cases was kept short in order to reduce the implementation burden in ATP; with the proposed method, the length and discretization of the line did not represent a problem.

Highly Nonuniform Transmission Line
A highly nonuniform three-phase transmission line was considered. The nonuniformity was included by the sagging of the conductors as shown in Figure 6.     Additionally, a comparison at the middle of the line between the transient voltage waveforms obtained with the proposed method and results from ATP simulations is presented (Figure 8). The simulations performed with ATP required a time step 10 times smaller than the proposed method to achieve similar results.  This example was simulated with the proposed method considering 20 and 40 chain matrices. Figure 9 presents a comparison of voltage measurements at the middle point of the transmission line from simulations considering 20, 40, and 80 chain matrices. From this comparison, it can be observed that the curves from the three simulations, in general, have the same shape (curves for 40 and 80 chain matrices are overlapped). However, the curve corresponding to the simulation performed with 20 chain matrices presents some oscillations in comparison with the other two. As it can be seen with the curves corresponding to 40 and 80 chain matrices, with an increase in the number of chain matrices used in the simulations, the curves become smoother; however, a further increase in the number of chain matrices used is barely noticeable. This indicates that although the method's accuracy is dependent on the number of chain matrices considered, it does not require a large number of chain matrices to achieve good results.  This example was simulated with the proposed method considering 20 and 40 chain matrices. Figure 9 presents a comparison of voltage measurements at the middle point of the transmission line from simulations considering 20, 40, and 80 chain matrices. From this comparison, it can be observed that the curves from the three simulations, in general, have the same shape (curves for 40 and 80 chain matrices are overlapped). However, the curve corresponding to the simulation performed with 20 chain matrices presents some oscillations in comparison with the other two. As it can be seen with the curves corresponding to 40 and 80 chain matrices, with an increase in the number of chain matrices used in the simulations, the curves become smoother; however, a further increase in the number of chain matrices used is barely noticeable. This indicates that although the method's accuracy is dependent on the number of chain matrices considered, it does not require a large number of chain matrices to achieve good results. This example was simulated with the proposed method considering 20 and 40 chain matrices. Figure 9 presents a comparison of voltage measurements at the middle point of the transmission line from simulations considering 20, 40, and 80 chain matrices. From this comparison, it can be observed that the curves from the three simulations, in general, have the same shape (curves for 40 and 80 chain matrices are overlapped). However, the curve corresponding to the simulation performed with 20 chain matrices presents some oscillations in comparison with the other two. As it can be seen with the curves corresponding to 40 and 80 chain matrices, with an increase in the number of chain matrices used in the simulations, the curves become smoother; however, a further increase in the number of chain matrices used is barely noticeable. This indicates that although the method's accuracy is dependent on the number of chain matrices considered, it does not require a large number of chain matrices to achieve good results.

Sequential Energization
A 5 km transmission line was used in this case. The nonuniformity was introduced by means of the line's sag: there was a transmission tower every 500 m, the conductors' height was maximum at the towers (20 m), and the height was minimum at the midspan between towers (15 m). The consideration for the height's variation as well as the horizontal separation between conductors was the same as in the previous example. The line was connected to an AC voltage source (1 p.u.) at the sending node through a three-phase switch. The switch's poles operated in a sequential manner with closing times of 3, 6, and 9 ms in an ABC sequence. The line's receiving node was left open.
The transient voltage profiles of phases A and B are presented in Figures 10 and 11, respectively. The comparison with ATP simulations is shown in Figure 12.

Sequential Energization
A 5 km transmission line was used in this case. The nonuniformity was introduced by means of the line's sag: there was a transmission tower every 500 m, the conductors' height was maximum at the towers (20 m), and the height was minimum at the midspan between towers (15 m). The consideration for the height's variation as well as the horizontal separation between conductors was the same as in the previous example. The line was connected to an AC voltage source (1 p.u.) at the sending node through a three-phase switch. The switch's poles operated in a sequential manner with closing times of 3, 6, and 9 ms in an ABC sequence. The line's receiving node was left open.
The transient voltage profiles of phases A and B are presented in Figures 10 and 11, respectively. The comparison with ATP simulations is shown in Figure 12.

Sequential Energization
A 5 km transmission line was used in this case. The nonuniformity was introduced by means of the line's sag: there was a transmission tower every 500 m, the conductors' height was maximum at the towers (20 m), and the height was minimum at the midspan between towers (15 m). The consideration for the height's variation as well as the horizontal separation between conductors was the same as in the previous example. The line was connected to an AC voltage source (1 p.u.) at the sending node through a three-phase switch. The switch's poles operated in a sequential manner with closing times of 3, 6, and 9 ms in an ABC sequence. The line's receiving node was left open.
The transient voltage profiles of phases A and B are presented in Figures 10 and 11, respectively. The comparison with ATP simulations is shown in Figure 12.

Surge Arrester Operation
This example presents the operation of surge arresters during a direct lightning strike; due to their v-i characteristic, the arresters were modeled as nonlinear elements as described in Section 2.3. The same line configuration presented in Section 3.1 was considered. The line was excited by a direct lightning strike (1.2/50 µs), the impact point was at the phase A of the sending node, and the line was impedance-matched at both ends in order to avoid reflections. The injected current to the line was approximated by a double exponential current source defined as: where τ1 = 68.199 µs, τ2 = 0.405 µs, and I0 = 10.37 kA. Surge arresters were connected to the three phases at the receiving node. The nonlinear characteristic of the arresters was simulated by means of a five-segment piecewise-linear approximation. Table 1 presents the v-i coordinates of this representation.

Surge Arrester Operation
This example presents the operation of surge arresters during a direct lightning strike; due to their v-i characteristic, the arresters were modeled as nonlinear elements as described in Section 2.3. The same line configuration presented in Section 3.1 was considered. The line was excited by a direct lightning strike (1.2/50 µs), the impact point was at the phase A of the sending node, and the line was impedance-matched at both ends in order to avoid reflections. The injected current to the line was approximated by a double exponential current source defined as: where τ1 = 68.199 µs, τ2 = 0.405 µs, and I0 = 10.37 kA. Surge arresters were connected to the three phases at the receiving node. The nonlinear characteristic of the arresters was simulated by means of a five-segment piecewise-linear approximation. Table 1 presents the v-i coordinates of this representation.

Surge Arrester Operation
This example presents the operation of surge arresters during a direct lightning strike; due to their v-i characteristic, the arresters were modeled as nonlinear elements as described in Section 2.3. The same line configuration presented in Section 3.1 was considered. The line was excited by a direct lightning strike (1.2/50 µs), the impact point was at the phase A of the sending node, and the line was impedance-matched at both ends in order to avoid reflections. The injected current to the line was approximated by a double exponential current source defined as: where τ 1 = 68.199 µs, τ 2 = 0.405 µs, and I 0 = 10.37 kA.
Surge arresters were connected to the three phases at the receiving node. The nonlinear characteristic of the arresters was simulated by means of a five-segment piecewise-linear approximation. Table 1 presents the v-i coordinates of this representation. First, a simulation was performed without surge arresters connected to the line. Figure 13 presents the transient voltage profile along phase A computed in this simulation. On the other hand, Figure 14 shows the voltage profile along phase A when the arresters were connected at the receiving node of the line. By comparing Figures 13 and 14, the influence that the surge arresters' operation had on the magnitude of the voltages along the line is easily observed. The voltage along the line in Figure 14 was considerably lower in comparison with the transient profile in Figure 13. Additionally, a comparison is presented with results obtained from ATP simulations (Figure 15). In a similar way to the example in Section 3.1, the time step required by ATP was 10 times smaller than the proposed method in order to obtain similar results.  First, a simulation was performed without surge arresters connected to the line. Figure 13 presents the transient voltage profile along phase A computed in this simulation. On the other hand, Figure 14 shows the voltage profile along phase A when the arresters were connected at the receiving node of the line. By comparing Figures 13 and 14, the influence that the surge arresters' operation had on the magnitude of the voltages along the line is easily observed. The voltage along the line in Figure 14 was considerably lower in comparison with the transient profile in Figure 13. Additionally, a comparison is presented with results obtained from ATP simulations (Figure 15). In a similar way to the example in Section 3.1, the time step required by ATP was 10 times smaller than the proposed method in order to obtain similar results.    First, a simulation was performed without surge arresters connected to the line. Figure 13 presents the transient voltage profile along phase A computed in this simulation. On the other hand, Figure 14 shows the voltage profile along phase A when the arresters were connected at the receiving node of the line. By comparing Figures 13 and 14, the influence that the surge arresters' operation had on the magnitude of the voltages along the line is easily observed. The voltage along the line in Figure 14 was considerably lower in comparison with the transient profile in Figure 13. Additionally, a comparison is presented with results obtained from ATP simulations (Figure 15). In a similar way to the example in Section 3.1, the time step required by ATP was 10 times smaller than the proposed method in order to obtain similar results.

Discussion
Transient simulation of nonuniform transmission lines using traditional software (such as ATP in this paper) is a challenging task. This difficulty is due to the fact that the nonuniformities are typically approximated by cascade-connecting several small line segments. Such approximation requires the use of a very small time step in the simulation (at least five times smaller than the traveling time of the line segments [6]), which may translate into excessive simulation times and the saturation of the software's available memory, as was the case in this work. Additionally, the connection of several line segments can introduce errors in the simulation results [2].
In contrast, the proposed method is able to accurately compute transient profiles along nonuniform transmission lines with substantially larger time steps in comparison to ATP, resulting in a better memory usage as expected from the findings reported in [16] with regards to the use of the INLT algorithm. As it can be observed in Figures 7, 10, 11 and 14, the main advantage of the proposed method is the fact that it allows visualizing the voltage and current transient behavior along the line and not only at its ends, which can be advantageous when designing lines with a high level of nonuniformities and cannot be easily done using traditional simulation software.
The precision of the proposed method was validated by comparisons with the results obtained with ATP simulations. In general, there was a very good level of agreement between the results from both methods, as it can be seen by the mean relative difference presented in Tables 2-4. There is a slight difference in the comparison of results presented in Figure 12 corresponding to the simulation of sequential pole closure. This is attributed to the previously mentioned ATP limitations and the difference in the time step used in each method, which can result in deviations in the operation time of the switch model, leading to the observed variations.

Discussion
Transient simulation of nonuniform transmission lines using traditional software (such as ATP in this paper) is a challenging task. This difficulty is due to the fact that the nonuniformities are typically approximated by cascade-connecting several small line segments. Such approximation requires the use of a very small time step in the simulation (at least five times smaller than the traveling time of the line segments [6]), which may translate into excessive simulation times and the saturation of the software's available memory, as was the case in this work. Additionally, the connection of several line segments can introduce errors in the simulation results [2].
In contrast, the proposed method is able to accurately compute transient profiles along nonuniform transmission lines with substantially larger time steps in comparison to ATP, resulting in a better memory usage as expected from the findings reported in [16] with regards to the use of the INLT algorithm. As it can be observed in Figures 7, 10, 11 and 14, the main advantage of the proposed method is the fact that it allows visualizing the voltage and current transient behavior along the line and not only at its ends, which can be advantageous when designing lines with a high level of nonuniformities and cannot be easily done using traditional simulation software.
The precision of the proposed method was validated by comparisons with the results obtained with ATP simulations. In general, there was a very good level of agreement between the results from both methods, as it can be seen by the mean relative difference presented in Tables 2-4. There is a slight difference in the comparison of results presented in Figure 12 corresponding to the simulation of sequential pole closure. This is attributed to the previously mentioned ATP limitations and the difference in the time step used in each method, which can result in deviations in the operation time of the switch model, leading to the observed variations.

Conclusions
This paper describes a frequency domain method to compute transient voltage and current profiles along nonuniform multiconductor transmission lines, where the nonuniformities along the line are introduced in the model by means of the cascaded connection of chain matrices. The method can incorporate nonlinear and time-dependent elements by using the superposition principle. The profiles obtained provide useful information to locate possible overvoltages at interior points along the transmission line, in contrast to traditional methods that only provide information at the line's ends. This information can be used as a helpful instrument in the insulation coordination design of transmission lines, as well as an educational tool in electrical engineering courses.
The results computed with the described method were compared with time-domain simulations using the well-known software ATP. In all of the comparisons, a high level of agreement was observed, demonstrating that the proposed method has a high level of accuracy.
It is worth mentioning that, although the results from both methods were very similar, the ATP simulations required substantially more time samples (10 to 50 times more samples) to achieve such results, which can result in a large computational burden. Additionally, in order to obtain measurements at interior points of the transmission line and to include nonuniformities along the transmission line in the ATP simulations, it is necessary to subdivide the line into many line segments, which can be a time-consuming process and can lead to error accumulation. Because of these issues, the proposed frequency domain method is a superior alternative to analyze the transient behavior of a transmission line when the voltage at interior points of the line is of interest.
The relevance of the proposed method lies in its strong potential for application in the accurate fault location and insulation coordination of transmission systems with prominent nonuniformities, which are commonly encountered in many power systems around the world.