**1. 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–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.

#### **2. 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.

#### *2.1. 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:

$$
\begin{bmatrix} V\_L(s) \\ I\_L(s) \end{bmatrix} = \begin{bmatrix} A & B \\ \mathbf{C} & D \end{bmatrix} \begin{bmatrix} V\_0(s) \\ I\_0(s) \end{bmatrix} \tag{1}
$$

where

$$\begin{aligned} A &= \cosh(\mathbf{y}L) \\ B &= -\mathbf{Z}\_0 \sinh(\mathbf{y}L) \\ \mathbf{C} &= \mathbf{Y}\_0 \sinh(\mathbf{y}L) \\ D &= -\cosh(\mathbf{y}L) \\ \mathbf{y} &= \sqrt{\mathbf{Z}\mathbf{Y}} \end{aligned} \tag{2}$$

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*<sup>ω</sup>, where ω is given by <sup>2</sup>π*f*.

**Figure 1.** Uniform transmission line representation. *V*0 and *VL* are the voltages of the sending and receiving node, and *I*0 and *IL* are the injected currents at the sending and receiving node, respectively, and *L* is the line's length.

By changing the direction of *IL* in Figure 1, Equation (1) is modified in the following manner:

$$
\begin{bmatrix} V\_L(s) \\ I\_L(s) \end{bmatrix} = \begin{bmatrix} A & B \\ -\mathbf{C} & -\mathbf{D} \end{bmatrix} \begin{bmatrix} V\_0(s) \\ I\_0(s) \end{bmatrix} \tag{3}
$$

or in a compact form:

$$
\begin{bmatrix} V\_L(s) \\ I\_L(s) \end{bmatrix} = \Phi \begin{bmatrix} V\_0(s) \\ I\_0(s) \end{bmatrix} \tag{4}
$$

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.** Cascaded connection of chain matrices, where **Φ***<sup>n</sup>*, *V<sup>n</sup>*, and *In* are the chain matrix, sending voltage, and sending current of the *n*-th cascade-connected line, respectively.

Figure 2 can also be interpreted as one transmission line subdivided into several smaller line segments, where each segmen<sup>t</sup> 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 segmen<sup>t</sup> and the chain matrix of the same segmen<sup>t</sup> can be used to compute the voltage and current of the next segment:

$$
\begin{bmatrix} V\_1(s) \\ I\_1(s) \end{bmatrix} = \Phi\_1 \begin{bmatrix} V\_0(s) \\ I\_0(s) \end{bmatrix}
$$

$$
\begin{bmatrix} V\_2(s) \\ I\_2(s) \end{bmatrix} = \Phi\_2 \begin{bmatrix} V\_1(s) \\ I\_1(s) \end{bmatrix} = \begin{bmatrix} \Phi\_2 \Phi\_1 \end{bmatrix} \begin{bmatrix} V\_0(s) \\ I\_0(s) \end{bmatrix} \tag{5}
$$

or in a general way:

$$\begin{bmatrix} V\_N(s) \\ I\_N(s) \end{bmatrix} = \Phi\_N \Phi\_{N-1} \dots \Phi\_3 \Phi\_2 \Phi\_1 \begin{bmatrix} V\_0(s) \\ I\_0(s) \end{bmatrix} \tag{6}$$

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 segmen<sup>t</sup> 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:

$$
\Phi\_{\rm FL} = \Phi\_N \Phi\_{N-1} \dots \Phi\_3 \Phi\_2 \Phi\_1 = \begin{bmatrix}
\Phi\_{\rm FL11} & \Phi\_{\rm FL12} \\
\Phi\_{\rm FL21} & \Phi\_{\rm FL22}
\end{bmatrix} \tag{7}
$$

$$
\begin{bmatrix} I\_0(s) \\ I\_L(s) \end{bmatrix} = \begin{bmatrix} \Upsilon\_{\mathfrak{s}\mathfrak{s}} & -\Upsilon\_{\mathfrak{s}\mathfrak{r}} \\ -\Upsilon\_{\mathfrak{s}\mathfrak{s}} & \Upsilon\_{\mathfrak{r}\mathfrak{r}} \end{bmatrix} \begin{bmatrix} V\_0(s) \\ V\_L(s) \end{bmatrix} \tag{8}
$$

where

$$\begin{array}{c} \mathbf{Y}\_{ss} = -\boldsymbol{\Phi}\_{\operatorname{FL12}}^{-1} \boldsymbol{\Phi}\_{\operatorname{FL11}}\\ \mathbf{Y}\_{\mathcal{S}\mathcal{T}} = -\boldsymbol{\Phi}\_{\operatorname{FL12}}^{-1} \\ \mathbf{Y}\_{\mathcal{S}\mathcal{T}} = \boldsymbol{\Phi}\_{\operatorname{FL21}} - \boldsymbol{\Phi}\_{\operatorname{FL22}} \boldsymbol{\Phi}\_{\operatorname{FL12}}^{-1} \boldsymbol{\Phi}\_{\operatorname{FL11}}\\ \mathbf{Y}\_{\mathcal{T}\mathcal{T}} = -\boldsymbol{\Phi}\_{\operatorname{FL22}} \boldsymbol{\Phi}\_{\operatorname{FL12}}^{-1} \end{array} \tag{9}$$

*V*0(*s*) is obtained by solving (8) for the voltages vector. *I*0(*s*) is computed as follows:

$$I\_0(s) = \mathcal{Y}\_{ss} V\_0(s) - \mathcal{Y}\_{rr} V\_L(s) \tag{10}$$

#### *2.2. Modeling of Time-Varying Elements*

It can be di fficult 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 e fficient 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:

$$\mathbf{V}\_{\rm NL}(\mathbf{s}) = \begin{bmatrix} \ \mathbf{V}\_{\rm E}(\mathbf{s}) & \ \mathbf{V}\_{\rm 0}(\mathbf{s}) & \ \mathbf{V}\_{\rm L}(\mathbf{s}) \ \end{bmatrix}^{\rm T} \tag{11}$$

and can be obtained from the following expression [11]:

$$\mathbf{V}\_{\rm NL}(s) = \mathbf{Y} \mathfrak{b} \mathfrak{u} \mathfrak{s}\_0^{-1} \mathbf{I}\_{\rm N0} + \mathbf{Y} \mathfrak{b} \mathfrak{u} \mathfrak{s}\_1^{-1} \mathbf{I}\_{\rm N1} \tag{12}$$

**Figure 3.** (**a**) Transmission line before the closing operation of a switch connected at the sending node. (**b**) Transmission line after the closing operation of a switch connected at the sending node.

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 *IN*1 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:

$$\begin{array}{c|ccccc} & Y\_{11} & -Y\_{12} & -Y\_{13} \\ & -Y\_{21} & Y\_{22} & -Y\_{23} \\ & -Y\_{31} & -Y\_{32} & Y\_{33} \\ \end{array} \tag{13}$$

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: *Y*12 = *Y*21 = *Yswitcho*, where *Yswitcho* is the open switch's admittance between its poles, ideally *Yswitcho* = 0. *Y*11 = *Ys* + *Yswtcho*, with *Ys* being the admittance connected at the sending node, and *Y*22 = *YLL* + *Yswitcho* where *YLL* represents the self-admittance of the transmission line. The matrix *Ybus*1 (Figure 3b) can be built in a similar way, but replacing *Yswitcho* by *Yswitchc*, 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:

$$\mathbf{V}\_{\rm NL}(\mathbf{s}) = \mathbf{Y} \mathbf{b} \mathbf{u} \mathbf{s}\_0^{-1} \mathbf{I}\_{\rm N0} + \sum\_{i=1}^{m} \mathbf{Y} \mathbf{b} \mathbf{u} \mathbf{s}\_i^{-1} \mathbf{I}\_{\rm Ni} \tag{14}$$

where *Ybusi* and *IN*i 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 *VL*(*s*) from (11).

#### *2.3. 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.

**Figure 4.** Five-segment piece-wise approximation of a nonlinear *v-i* characteristic curve [11].

The element *Rn* of the *n*-th linear element from Figure 4 and the voltage *Vn* are computed as follows:

$$R\_n = \frac{V\_{refn+1} - V\_{refn}}{I\_{refn+1} - I\_{refn}} \tag{15}$$

$$V\_n = \left(-\mathcal{R}\_{n+1} I\_{refn}\right) + V\_{refn+1} \tag{16}$$

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.

**Figure 5.** Circuital representation of the piecewise approximation for the characteristic curve of a nonlinear element presented in Figure 5 [11].

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 *Rn*, that is, the slope of the *n*-th linear segmen<sup>t</sup> in the approximation of Figure 4:

$$R\_{\rm H} = \frac{R\_{\rm n-1} R\_{\rm xm}}{R\_{\rm n-1} + R\_{\rm xm}} \tag{17}$$

Finally, by solving (17) for *Rxn*, the value of such resistance can be computed:

$$R\_{\mathbf{x}\_n} = \frac{R\_{n-1}R\_n}{R\_{n-1} - R\_n} \tag{18}$$

*Vxn* is computed as:

$$V\_{\mathbf{x}\_n} = \frac{R\_{n-1}V\_n - V\_{n-1}R\_n}{R\_{n-1} - R\_n} \tag{19}$$

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.

#### *2.4. 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:

$$f(t) = \text{Re}\left\{\frac{e^{\epsilon t}}{2\pi} \int\_0^\infty F(s)e^{j\omega t} d\alpha\right\} \tag{20}$$

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:

$$f\_n = f(n\Delta t)F\_m = F(c + j(2m+1)\Delta \omega) \tag{21}$$

where *n*, *m* = 0, 1, 2, ... , *N* − 1 and *c*hat reduces the aliasing errors of the algorithm; in this work, *c* is defined as *c* = 2Δ<sup>ω</sup>.

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:

$$T = \frac{\pi}{\Delta w} \tag{22}$$

and the following relations can be established:

$$\begin{aligned} \Delta t &= \frac{T}{N} \\ \Delta \omega &= \frac{\Omega}{2N} = \frac{\pi}{T} \end{aligned} \tag{23}$$

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:

$$f\_n = \operatorname{Re}\left\{ \mathbf{C}\_n \left[ \sum\_{m=0}^{N-1} F\_m \sigma\_m \exp\left(\frac{j2\pi m}{N}\right) \right] \right\} \\ \text{for } n, m = 0, 1, 2, \dots, N - 1 \tag{24}$$

where:

$$C\_{n} = \frac{2\Delta\omega}{\pi} \exp\left(cn\Delta t + \frac{j\pi n}{N}\right) \tag{25}$$

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.

#### **3. 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.

#### *3.1. 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. The line was excited by an ideal unit step voltage source (1 p.u.) connected at phase A of the sending node, while all of the other nodes were left open.

**Figure 6.** Three-phase non-uniform transmission line used in the example in Section 3.1. A second-degree polynomial equation was utilized to approximate the conductors' height, as presented in [23]. There was a 10 m separation between the conductors.

Figure 7 presents the transient voltage profile along phase A computed with the proposed method; in order to include the nonuniformities along the transmission line, it was subdivided into 80 chain matrices. This figure illustrates the propagation of the traveling waves along the line and how these waves were reflected when they reached the receiving end. It can also be observed in the transient profiles that there were periods where negative voltages appeared at some points along the line, which is not a commonly observed phenomenon when similar simulations with uniform lines are performed. 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.

**Figure 7.** Computed voltage profile along phase A (example 3.1).

**Figure 8.** Comparison of transient voltage waveforms at two specific distances between the computed results with the proposed method (solid lines) and results obtained with ATP simulations (dashed lines). The solid blue and dashed red lines correspond to voltage measurements at 150 m from the sending node; the solid green and the dashed orange lines were obtained at 450 m. The oscillations observed in the results from the ATP simulations are attributed to the error accumulation due to the discretization of the transmission line [2].

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.

**Figure 9.** Comparison of transient voltage waveforms at the middle point of the line computed with different numbers of chain matrices. The curve computed with 20 chain matrices is not as smooth as the ones computed with 40 and 80.

## *3.2. 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.

**Figure 10.** Voltage profile computed along phase A (example 3.2).

**Figure 11.** Transient voltage profile along phase B in example 3.2. The complete voltage profile induced in phase B due to the energization of phase A can be observed. This type of detailed information cannot be easily achieved using ATP-like simulation software.

**Figure 12.** Comparison between the transient voltage computed with the proposed method (blue and green lines) and those obtained with ATP simulations (red and orange lines) in example 3.2. The comparisons are made at the middle of the transmission line (2500 m from the sending node) in phases A and B. A good level of agreemen<sup>t</sup> can be observed, but the ATP simulations needed a time step 50 times smaller than the proposed method in order to achieve such results.

#### *3.3. 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:

$$i(t) = I\_0 \left( e^{-\frac{t}{\tau\_1}} - e^{-\frac{t}{\tau\_2}} \right) \tag{26}$$

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.


**Table 1.** Nonlinear *v-i* characteristic of the arresters [27].

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.

**Figure 13.** Computed transient voltage profile along phase A when there were no surge arresters connected at the receiving end of the line (example 3.3).

**Figure 14.** Computed transient voltage profile along phase A with surge arresters connected at the receiving node (example 3.3). The voltage magnitude significantly decreased along the transmission line compared to the transient profile presented in Figure 13.

**Figure 15.** Voltage comparison between the results from the proposed method (solid lines) and those obtained from ATP simulations (dashed lines) for example 3.3. The comparisons were made at the middle point of the line (300 m from the sending node) for phases A, B, and C.
