Analysis of Metal Oxide Varistor Arresters for Protection of Multiconductor Transmission Lines Using Unconditionally-Stable Crank–Nicolson FDTD

Abstract

Surge arresters may represent an efficient choice for limiting lightning surge effects, significantly reducing the outage rate of power lines. The present work firstly presents an efficient numerical approach suitable for insulation coordination studies based on an implicit Crank–Nicolson finite difference time domain method; then, the IEEE recommended surge arrester model is reviewed and implemented by means of a local implicit scheme, based on a set of non-linear equations, that are recast in a suitable form for efficient solution. The model is proven to ensure robustness and second-order accuracy. The implementation of the arrester model in the implicit Crank–Nicolson scheme represents the added value brought by the present study. Indeed, its preserved stability for larger time steps allows reducing running time by more than $60\%$ compared to the well-known finite difference time domain method based on the explicit leap-frog scheme. The reduced computation time allows faster repeated solutions, which need to be looked for on assessing the lightning performance (randomly changing, parameters such as peak current, rise time, tail time, location of the vertical leader channel, phase conductor voltages, footing resistance, insulator strength, etc. would need to be changed thousands of times).

1. Introduction

Overhead power lines directly affect the reliable and safe operation of power grids. In consideration of their long distances and large capacity power transmission, insulation coordination of Multiconductor Transmission Lines (MTLs) needs to include protection for fast and slow front transient overvoltage stresses caused by lightning strikes [1,2,3,4] or switching since they can be critical for insulation and pose a danger to connected equipments [5]. Surge protection is a main concern in the design of wind farms [6,7] and high-voltage (HV) substations [8].

According to recent practices of insulation coordination based on a statistical approach of electrical stress and strength [5,9], direct/indirect lightning strikes can cause an increasing probability of insulator string flashovers when the voltage transient surges exceed the Insulation Level (IL). Protection practice usually relies on the use of overhead shield wires [10] that are grounded at the towers [11,12]. However, the grounding resistance must be kept as low as possible because the discharge of lightning current through the ground system results in a Ground Potential Rise (GPR) [13], which can cause the voltage across the insulator strings to be higher than their IL, so that back-flashover of the insulation will occur [14].

According to recent studies [15], overvoltage limiting devices consisting of surge arresters or Metal Oxide Varistor (MOV) [16] arresters can be effectively located along the power lines. Optimized surge arrester allocation based on genetic algorithm or other techniques have been proposed (e.g., [17]). In particular, MOVs may represent an efficient choice for limiting lightning and line-switching surges, significantly reducing the outage rate. It is well-documented [18,19] that the use of MOVs may offer additional advantages, which reduce the need for complex solutions to limit switching overvoltages, e.g., controlled closing of circuit breakers equipped with preinsertion resistors.

Nowadays [20,21], gapless lightweight metal-oxide type line surge arresters (LSA) represent a very promising solution since they can improve the protection level, offering special features inherent gapless structure, e.g., quick response to surges, high energy dissipation capability, and freedom from pollution, that may introduce spurious resistive components [22].

The metal-oxide granular structure and chemical properties confer to MOVs their well-known nonlinear V–I characteristic [23,24,25]. The varistor behavior allows the device to limit the level of the voltage, avoiding undesirable values, and to dissipate part of the energy associated with the travelling surge by draining current to the ground. However, if the energy absorbed by the MOVs exceeds their energy handling capability that is partially described by the line discharge class (event that can occur under direct strike), they will be damaged. For this reason, the authors of [26] concluded that a correct coordination among MOV withstand capability, number of shield wires, and grounding resistance of arresters is necessary to limit permanent damage to surge arresters and obtain an effective lightning protection.

Several models have been proposed to describe the MOV behavior [27], namely IEEE model recommended by Working Group (WG) 3.4.11 [28], Pinceti model [29], Valsalal model [30], and Fernandez model [31]. Different characteristics are associated with different construction features, the choice of which is linked to the required performance. Appropriate modeling of MOV dynamic characteristics is certainly of paramount importance for insulation coordination studies and arrester location [17].

Improving modeling techniques of protection devices and consideration of protection in analytical tools are activities still being researched. Commercial codes (e.g., EMTP-RV and ATP-EMTP) based on approximated line models (e.g., JMarti or Bergeron model) and numerical tools based on the Finite Difference Time Domain (FDTD) method are widely available, but their accuracy is still debated.

The present work presents a novel numerical methodology to study the dynamic behavior of MOVs in power lines. The analysis is carried out in the Time-Domain (TD) [32] applying an implicit Crank–Nicolson Finite Difference (CN-FD) method, which allows achieving reliable results with a significant reduction of simulation time. The IEEE recommended model is revised and efficiently implemented, for the first time, with a local implicit scheme by using a set of non-linear equations that ensures robustness and second-order accuracy. Comparison with classical explicit FDTD based on a leap-frog marching-on-time scheme [33] shows the accuracy of the proposed method and how it may run substantially faster if the Courant stability condition is relaxed. In addition, we propose a smart strategy to include the MOV model inside the overall implicit scheme.

Using the proposed methodology, the paper carries out an analysis on the performance of gapless lightweight LSA installed on power transmission lines. The methodology proves to be particularly well suited to provide guidance on best allocation of MOVs in terms of their number and position, and on the coordination between MOVs and grounding resistance of the tower earthing system. The technique is promising in insulation coordination studies and the authors believe that it can be also applied to distribution lines, where novel distribution network architectures are nowadays required.

The rest of the paper is organized as follows. In Section 2, we present the novel numerical methodology based on the implicit Crank–Nicolson Finite Difference method and explain how the IEEE MOV model can be efficiently inserted. In Section 3, we present the insulation coordination study of a 150 kV line and discuss the main features of the proposed numerical methodology. Finally, we draw some conclusions in Section 4.

2. Numerical Formulation

Under the quasi-transverse electromagnetic (TEM) transmission line approximation, the equations describing the propagation on an overhead multiconductor transmission line (MTL) in the presence of a finite, conductive ground and corona read (we use lowercase boldface letters for vector quantities and capital boldface letters for matrix quantities):

$$\begin{array}{cc}\hfill -\frac{\partial \mathbf{v}\left(x,t\right)}{\partial x}& ={\int }_0^t{\mathbf{Z}}_{\zeta }^{\prime }\left(t-\tau \right)\frac{\partial \mathbf{i}}{\partial \tau }d\tau +{\mathbf{L}}_{\mathrm{ext}}^{\prime }\frac{\partial \mathbf{i}\left(x,t\right)}{\partial t}\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill -\frac{\partial \mathbf{i}\left(x,t\right)}{\partial x}& ={\mathbf{C}}^{\prime }\frac{\partial \mathbf{v}\left(x,t\right)}{\partial t},\hfill \end{array}$$

(1b)

where $\mathbf{v}$ and $\mathbf{i}$ are the voltages and currents at coordinate x along the line and at time instant t, respectively; ${\mathbf{L}}_{\mathrm{ext}}^{\prime }$ is the external line inductance per unit length (p.u.l.); ${\mathbf{Z}}_{\zeta }^{\prime }$ is the transient impedance which is usually introduced in the TD to include the frequency-dependent behavior of ground (ground transient impedance [34,35]) and wires (internal transient impedance [36]); and ${\mathbf{C}}^{\prime }$ is the p.u.l. geometric capacitance. The matrix ${\mathbf{C}}^{\prime }$ can be replaced with a non-linear time-dependent dynamic capacitance ${\mathbf{C}}_{\mathrm{dyn}}$ to take into account for corona effect [37,38], which is neglected in the present paper. In Equation (1), the term corresponding to the so-called ground admittance is neglected. In fact, this is a reasonable approximation for typical overhead lines and for the frequency ranges of interest, namely below 10 MHz. Moreover, we neglect the effect of the p.u.l. conductance $\mathbf{G}$, since transversal losses are usually associated with the corona effect and accounted for through the dynamic capacitance ${\mathbf{C}}_{\mathrm{dyn}}$, as aforementioned. It should be observed that the inclusion of corona effect is particularly challenging [38] and several models are available. Nevertheless, from a numerical standpoint, the dynamic change of the p.u.l. time-dependent capacitance affects the propagation velocity of the surge along the line and could result in instabilities in an explicit method if the Courant stability condition is no more fulfilled. This problem is not encountered in an implicit unconditionally stable method, where only a careful control of the space sampling is requested.

The overhead line is assumed to be made of $N_{\mathrm{w}}$ wires (see Figure 1), i.e., $N_{\mathrm{p}}$ phase conductors and $N_{\mathrm{s}}$ overhead shield wires (OHSWs), with $N_{\mathrm{w}}=N_{\mathrm{p}}+N_{\mathrm{s}}$. The length of the line is $L=N_{\mathrm{span}}\ell$, where ℓ is the pole-to-pole span and $N_{\mathrm{span}}$ is the number of sections. The OHSWs are grounded at some poles by a constant equivalent ground admittance $G_{\mathrm{g}}$. Propagation along the ground lead and its electrical parameters are neglected here [39].

From a practical point of view, we account for the effect of the sag $s_i$ of the ith conductor through an equivalent height $h_{\mathrm{e},i}$ of the wire computed as $h_{\mathrm{e},i}=h_i-\frac{2}{3}{s}_i$, where $h_i$ is the height of the ith conductor at the pole. Indeed, this approximation is valid for the frequencies of interest (below 10 MHz) and allows us to simplify the analysis, using the standard uniform transmission-line theory [40,41].

MOVs are assumed to be installed on one or more phases along the transmission line with the purpose of limiting both lightning and switching overvoltages under a specified IL. The MOVs should be chosen to let the residual voltages be lower than the dielectric strength of insulators, even for the highest strike currents. Nonetheless, laboratory test data of MOV discharges indicate that arresters have dynamic characteristics that play a fundamental role for studies involving fast front surges [23]. The induction of steep-fronted and high-magnitude transients can be enhanced by the discharge of MOVs, as well as through the common grounding system, and lead to severe failures of insulation of the non-protected phases whenever these devices are installed on just one or two phases. For these reasons, the frequency-dependent IEEE model recommended by WG 3.4.11 is considered here. It gives satisfactory results for discharge currents with a wide range of rise times, from 0.5 to 45 $\mathsf{\mu }$s. The main deadlocks when dealing with this model are its strongly nonlinear behavior, especially under fast transients, and the need for a careful and accurate selection of its constitutive parameters.

2.1. Crank–Nicolson Scheme

The Crank–Nicolson (CN) scheme [42,43,44,45] is second-order accurate in both space and time and unconditionally stable, being the size of the time step governed only by the accuracy and not by stability. The method is also equivalent to a single average alternating direction implicit (ADI) scheme [46].

The idea underlying the CN scheme is to apply centered differences in space and time, combined with an average in time [47]. To do so, we introduce a staggered mesh in space with constant discretization parameter $\Delta x$, where we seek voltages ${\mathbf{v}}_k$ at $N_{\Delta \mathrm{x}}+1$ integer mesh points $x_k$ (the mesh boundaries are voltage nodes) and currents ${\mathbf{i}}_{k+\frac{1}{2}}$ at $N_{\Delta \mathrm{x}}$ half-integer mesh points $x_{k+\frac{1}{2}}$ between the voltage points (we use the typical notation for finite differences as in [33]). Sampling both voltages and currents at the same time points $t^n$, the CN scheme starts by splitting each time step $\Delta t$ into two sub-steps $\Delta t/2$. In the first sub-step the explicit Forward Euler method is used to deduce the unknowns at time $\left(n+\frac{1}{2}\right)\Delta t$:

$$\begin{array}{cc}\hfill -\frac{{\mathbf{v}}_{k+1}^n-{\mathbf{v}}_k^n}{\Delta x}& ={\mathbf{L}}_{\mathrm{ext}}^{\prime }\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+\frac{1}{2}}-{\mathbf{i}}_{k+\frac{1}{2}}^n}{\Delta t/2}\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill -\frac{{\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k-\frac{1}{2}}^n}{\Delta x}& ={\mathbf{C}}^{\prime }\frac{{\mathbf{v}}_k^{n+\frac{1}{2}}-{\mathbf{v}}_k^n}{\Delta t/2}.\hfill \end{array}$$

(2b)

In the second sub-step, we use the implicit Backward Euler method to compute the unknowns at time $\left(n+1\right)\Delta t$:

$$\begin{array}{cc}\hfill -\frac{{\mathbf{v}}_{k+1}^{n+1}-{\mathbf{v}}_k^{n+1}}{\Delta x}& ={\mathbf{L}}_{\mathrm{ext}}^{\prime }\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t/2}\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill -\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k-\frac{1}{2}}^{n+1}}{\Delta x}& ={\mathbf{C}}^{\prime }\frac{{\mathbf{v}}_k^{n+1}-{\mathbf{v}}_k^{n+\frac{1}{2}}}{\Delta t/2}.\hfill \end{array}$$

(3b)

To our purpose, it may be more convenient to rearrange the two sub-steps in Equations (2) and (3) into a single time-step as

$$\begin{array}{cc}\hfill {\mathbf{i}}_{k+\frac{1}{2}}^{n+1}& ={\mathbf{i}}_{k+\frac{1}{2}}^n-\frac{\Delta t}{2}{\left[{\mathbf{L}}_{\mathrm{ext}}^{\prime }\right]}^{-1}\left(\frac{{\mathbf{v}}_{k+1}^{n+1}-{\mathbf{v}}_k^{n+1}}{\Delta x}+\frac{{\mathbf{v}}_{k+1}^n-{\mathbf{v}}_k^n}{\Delta x}\right)\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill {\mathbf{v}}_k^{n+1}& ={\mathbf{v}}_k^n-\frac{\Delta t}{2}{\left[{\mathbf{C}}^{\prime }\right]}^{-1}\left(\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k-\frac{1}{2}}^{n+1}}{\Delta x}+\frac{{\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k-\frac{1}{2}}^n}{\Delta x}\right).\hfill \end{array}$$

(4b)

The last two equations (which are valid for lossless MTLs) make evident that the CN scheme is based on a central difference on time $\left(n+\frac{1}{2}\right)\Delta t$ and on an interpolation of the space derivative term. This characteristic ensures second-order accuracy in both space and time.

This algorithm is considerably more intricate than the classical leap-frog FDTD based on the Yee scheme [48,49], since Equation (4) must be solved simultaneously, requiring the solution of a system of equations for each time step, which could be a costly procedure. The difficulty to handle this computational cost has been always a relevant obstacle to using the CN scheme in practice. Nonetheless, several considerations apply. Through suitable manipulations, Equations (4a) and (4b) may be combined to obtain:

$$\begin{array}{cc}\hfill -{\mathbf{v}}_k^{n+1}+\frac{\Delta x{\mathbf{L}}_{\mathrm{ext}}^{\prime }}{\frac{\Delta t}{2}}{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}+{\mathbf{v}}_{k+1}^{n+1}& =\frac{\Delta x{\mathbf{L}}_{\mathrm{ext}}^{\prime }}{\frac{\Delta t}{2}}{\mathbf{i}}_{k+\frac{1}{2}}^n-\left({\mathbf{v}}_{k+1}^n-{\mathbf{v}}_k^n\right)\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill -{\mathbf{i}}_{k-\frac{1}{2}}^{n+1}+\frac{\Delta x{\mathbf{C}}^{\prime }}{\frac{\Delta t}{2}}{\mathbf{v}}_k^{n+1}+{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}& =\frac{\Delta x{\mathbf{C}}^{\prime }}{\frac{\Delta t}{2}}{\mathbf{v}}_k^n-\left({\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k-\frac{1}{2}}^n\right);\hfill \end{array}$$

(5b)

We solve Equation (5) alternating a voltage node with a current node in the unknown column vector ${\mathbf{x}}^{n+1}$ and we obtain a linear system in the form:

$$\mathbf{A}{\mathbf{x}}^{n+1}={\mathbf{b}}^{n+1}$$

(6)

with

$$\begin{array}{cc}\hfill {\mathbf{x}}^{n+1}& =\left[{\mathbf{v}}_0^{n+1}{\mathbf{i}}_{\frac{1}{2}}^{n+1}{\mathbf{v}}_1^{n+1}\cdots {\mathbf{v}}_{N_{\Delta \mathrm{x}}}^{n+1}\right]\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill {\mathbf{b}}^{n+1}& =\left[{\mathbf{b}}_0^{\mathrm{V}}{\mathbf{b}}_{\frac{1}{2}}^{\mathrm{I}}{\mathbf{b}}_1^{\mathrm{V}}\cdots {\mathbf{b}}_{N_{\Delta \mathrm{x}}}^{\mathrm{V}}\right].\hfill \end{array}$$

(7b)

In Equation (6), $\mathbf{A}$ is a block tri-diagonal matrix (the size of each block is $N_{\mathrm{w}}\times N_{\mathrm{w}}$) whose band is $2N_{\mathrm{w}}-1$:

$$\mathbf{A}=\left[\begin{array}{ccccccc}{\mathbf{D}}_0^{\mathbf{V}}& \mathbf{I}& 0& 0& 0& \cdots & 0\\ -\mathbf{I}& {\mathbf{D}}_1^{\mathbf{I}}& \mathbf{I}& 0& 0& \cdots & 0\\ 0& -\mathbf{I}& {\mathbf{D}}_1^{\mathbf{V}}& \mathbf{I}& 0& \cdots & 0\\ ⋮& & & \ddots & & & ⋮\\ 0& & -\mathbf{I}& {\mathbf{D}}_k^{\mathbf{V}}& \mathbf{I}& & 0\\ ⋮& & & & & \ddots & ⋮\\ 0& \cdots & 0& 0& 0& -\mathbf{I}& {\mathbf{D}}_{N_{\Delta \mathrm{x}}}^{\mathbf{V}}\end{array}\right].$$

(8)

If the p.u.l. parameters are time-independent (i.e., we neglect corona effect that results in the time-dependent dynamic capacitance matrix ${\mathbf{C}}_{\mathrm{dyn}}\left(t\right)$), the matrix $\mathbf{A}$ is fixed and needs to be computed only once. Moreover, on each voltage/current node, two of the three blocks are identity (unitary) matrices $\mathbf{I}$. Consequently, the linear matrix system can be efficiently solved in band-storage form, using a parallelized or GPU CUDA solver whose performance scales with bandwidth (related to the number of wires $N_{\mathrm{w}}$) and not with matrix dimension [50].

The evaluation of the ground resistive-inductive effects follows from inserting into Equation (4a) the discrete counterpart of the ground transient impedance. To maintain a centered synchronized updating equation, we evaluate the loss term at time $\left(n+\frac{1}{2}\right)\Delta t$ as

$$\begin{array}{cc}-\left(\frac{{\mathbf{v}}_{k+1}^{n+1}-{\mathbf{v}}_k^{n+1}}{2\Delta x}+\frac{{\mathbf{v}}_{k+1}^n-{\mathbf{v}}_k^n}{2\Delta x}\right)& ={\mathbf{L}}_{\mathrm{ext}}^{\prime }\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k+\frac{1}{2}}^n}{\Delta t}+\hfill \\ & +\underset{0}{\overset{\left(n-\frac{1}{2}\right)\Delta t}{\int }}{\mathbf{Z}}_{\zeta }^{\prime }\left(t-\tau \right){\left.\frac{\partial \mathbf{i}}{\partial \tau }\right|}_{k+\frac{1}{2}}\mathrm{d}\tau +\underset{\left(n-\frac{1}{2}\right)\Delta t}{\overset{\left(n+\frac{1}{2}\right)\Delta t}{\int }}{\mathbf{Z}}_{\zeta }^{\prime }\left(t-\tau \right){\left.\frac{\partial \mathbf{i}}{\partial \tau }\right|}_{k+\frac{1}{2}}\mathrm{d}\tau .\hfill \end{array}$$

(9)

If we use a trapezoidal rule to approximate the second integral over the last time-step and we denote the first history integral as $\mathbf{H}(x,t)$, we end up with

$$\begin{array}{cc}-\left(\frac{{\mathbf{v}}_{k+1}^{n+1}-{\mathbf{v}}_k^{n+1}}{2\Delta x}+\frac{{\mathbf{v}}_{k+1}^n-{\mathbf{v}}_k^n}{2\Delta x}\right)& ={\mathbf{L}}_{\mathrm{ext}}^{\prime }{\displaystyle \frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k+\frac{1}{2}}^n}{\Delta t}}+\hfill \\ & +{\mathbf{H}}_{k+\frac{1}{2}}^{n-\frac{1}{2}}+{\displaystyle \frac{{\mathbf{Z}}_{\zeta }^{\prime }\left(0\right)}{2}}\left({\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k+\frac{1}{2}}^n\right)+{\displaystyle \frac{{\mathbf{Z}}_{\zeta }^{\prime }\left(\Delta t\right)}{2}}\left({\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k+\frac{1}{2}}^{n-1}\right).\hfill \end{array}$$

(10)

Hence, the presence of losses requires an additional term to be included in the second implicit sub-step, being the first explicit sub-step left unchanged:

$$\begin{array}{cc}\hfill -\frac{{\mathbf{v}}_{k+1}^{n+1}-{\mathbf{v}}_k^{n+1}}{\Delta x}& ={\mathbf{L}}_{\mathrm{ext}}^{\prime }\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t/2}+\hfill \\ \hfill & +2\left[{\mathbf{H}}_{k+\frac{1}{2}}^{n-\frac{1}{2}}+\frac{\Delta t}{2}\left({\mathbf{Z}}_{\zeta }^{\prime }\left(0\right)\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k+\frac{1}{2}}^n}{\Delta t}+{\mathbf{Z}}_{\zeta }^{\prime }\left(\Delta t\right)\frac{{\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k+\frac{1}{2}}^{n-1}}{\Delta t}\right)\right]\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill -\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k-\frac{1}{2}}^{n+1}}{\Delta x}& ={\mathbf{C}}^{\prime }\frac{{\mathbf{v}}_k^{n+1}-{\mathbf{v}}_k^{n+\frac{1}{2}}}{\Delta t/2}.\hfill \end{array}$$

(11b)

The second-order accuracy is still retained as well.

Equation (10) is an implicit equation which does not change the block tri-diagonal nature of the aforementioned scheme. Furthermore, it can be efficiently computed with a recursive technique based on Prony approximation of the transient impedance terms [32,38,51] as follows:

$$\begin{array}{cc}\hfill -{\mathbf{v}}_k^{n+1}+\frac{\Delta x{\mathbf{L}}_{\mathrm{ext}}^{\prime }}{\Delta t/2}{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}+{\mathbf{v}}_{k+1}^{n+1}& =\frac{\Delta x{\mathbf{L}}_{\mathrm{ext}}^{\prime }}{\Delta t/2}{\mathbf{i}}_{k+\frac{1}{2}}^n-\left({\mathbf{v}}_{k+1}^n-{\mathbf{v}}_k^n\right)+\hfill \\ \hfill & -2\Delta x\left[{\mathbf{H}}_{k\frac{1}{2}}^{n-\frac{1}{2}}+\left(-{\mathbf{Z}}_{\zeta }^{\prime }\left(0\right)\frac{{\mathbf{i}}_{k+\frac{1}{2}}^n}{2}+{\mathbf{Z}}_{\zeta }^{\prime }\left(\Delta t\right)\frac{{\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k+\frac{1}{2}}^{n-1}}{2}\right)\right]\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill -{\mathbf{i}}_{k-\frac{1}{2}}^{n+1}+\frac{\Delta x{\mathbf{C}}^{\prime }}{\Delta t/2}{\mathbf{v}}_k^{n+1}+{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}& =\frac{\Delta x{\mathbf{C}}^{\prime }}{\Delta t/2}{\mathbf{v}}_k^n-\left({\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k-\frac{1}{2}}^n\right).\hfill \end{array}$$

(12b)

2.2. MOV Numerical Model

The non-linear circuital model of surge arresters presented by WG 3.4.11 [28] is depicted in Figure 2. It consists of two R-L filters and two nonlinear resistors, one of which is in parallel with the capacitance C.

Under slow transients, the effect of the two R-L filters is negligible with respect to the resistance offered by the nonlinear branches; hence, $A_0$ and $A_1$ are essentially connected in parallel under the port voltage $v_{\mathrm{mov}}$. Under fast transients, the filters impedances cannot be neglected due to the wider range of frequencies excited, resulting in a significant voltage drop between the two nonlinear branches. The high-frequency current is forced by the $R_1$-$L_1$ filter to flow more in $A_0$ than in $A_1$. The requirement that the nonlinear resistor $A_0$ has higher voltage than $A_1$ at a given current, as shown in Figure 3, leads to the arrester model producing higher residual voltages at higher frequencies.

Several other models have been proposed in the literature [52]. The model shown in Figure 4a was proposed by Pinceti [29] who also developed a procedure for the parameters calculation. In [31], Fernandez et al. suggested another model which is also based on the IEEE model (see Figure 4b) and proposed an iterative trial and error procedure to compute the best model parameters. In [30], Valsalal et al. developed an arrester model for very fast transients (see Figure 4c) derived from a simplified IEEE model with the inclusion of the arrester block capacitance ($C_{\mathrm{bp}}$) and stray capacitance ($C_{\mathrm{g}}$).

A first approximation of the values to be given to the resistive, inductive, and capacitive parameters [28,53] is related to the arrester height $h_{sa}$ and to the number n of parallel columns of the MOV: $L_0=0.2h_{\mathrm{mov}}/n\left[\mathsf{\mu }\mathrm{H}\right]$, $R_0=100h_{\mathrm{mov}}/n\left[\Omega \right]$, $L_1=15h_{\mathrm{mov}}/n\left[\mathsf{\mu }\mathrm{H}\right]$, $R_1=65h_{\mathrm{mov}}/n\left[\Omega \right]$, $C=100n/h_{\mathrm{mov}}\left[\mathrm{pF}\right]$.

The voltage–current characteristics of the two nonlinear varistors $A_0$ and $A_1$ have been represented by several formulas or numerical models [28] determined from results of experimental investigations on arresters. These characteristics show a typical power-law trend $v=A_w{i}^k$, as shown in Figure 3, the voltage being expressed in per unit of the MOV maximum residual voltage when draining a 10 kA, 8/20 $\mathsf{\mu }$s impulse current. The parameters values must be optimized to reproduce the voltage–current characteristic of the MOV under analysis [54]. We used the MATLAB optimization toolbox to obtain the best parameters for our MOV model, i.e., those parameters that minimize the error between the transient voltage response of the model under 5 kA, 10 kA, 20 kA, 40 kA, 8/20 $\mathsf{\mu }$s impulse currents, and the measured voltage values in the data sheet.

Since the circuit has two nonlinear varistors, $A_0$ presenting higher voltage values than $A_1$ for a given current, as shown in Figure 3, the set of variables and, successively, of Kirchhoff’s equations must be properly chosen to ensure stability and robustness of the final discrete time-stepping scheme. The behavior of the circuit is captured by a set of equations that are formulated by combining the constitutive element equations and Kirchhoff’s Current and Voltage Laws (KCL and KVL). This results in a set of three simultaneous nonlinear first-order differential equations.

As concerns the choice of the variables, we choose three state variables: the two currents $i_{{\mathrm{L}}_0}$ and $i_{{\mathrm{L}}_1}$ through the inductors $L_0$ and $L_1$, respectively—that are mesh currents as well—and the voltage $e_0$ across the capacitor C and the varistor $A_0$. The voltage $e_0$ is used to control the varistor $A_0$, which is represented as a voltage-controlled current source (VCCS), and to find the current $i_{A0}$ through the nonlinear constitutive relation. Likewise, we represent the varistor $A_1$ as a VCCS: we derive the voltage $e_1$ through a KVL and we compute the current $i_{A1}$ through the nonlinear constitutive relation.

Starting with one KVL around Loop #1 and two KCLs through Cuts #1 and #2, as shown in Figure 2, and plugging the KVL around Loop #2 into the constitutive equation of $A_1$, we obtain:

$$\begin{array}{cc}\hfill & v_{\mathrm{mov}}-L_0\frac{\mathrm{d}{i}_{L_0}}{\mathrm{d}t}-e_0=0\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill & i_{L_0}+\frac{v_{\mathrm{mov}}-e_0}{R_0}-C\frac{\mathrm{d}{e}_0}{\mathrm{d}t}-i_{A0}\left({\mathrm{e}}_0\right)-i_{L_1}-\frac{L_1}{R_1}\frac{\mathrm{d}{i}_{L_1}}{\mathrm{d}t}=0\hfill \end{array}$$

(13b)

$$\begin{array}{cc}\hfill & i_{L_1}+\frac{L_1}{R_1}\frac{\mathrm{d}{i}_{L_1}}{\mathrm{d}t}-i_{A1}\left(e_0-L_1\frac{\mathrm{d}{i}_{L_1}}{\mathrm{d}t}\right)=0.\hfill \end{array}$$

(13c)

We observe that representing the varistors as VCCSs, the circuit is not degenerate.

Using a central difference combined with an average in time leads to the discrete version of Equation (13):

$$\begin{array}{cc}\hfill & L_0\frac{{\mathrm{i}}_{L_0}^{n+1}-{\mathrm{i}}_{L_0}^n}{\Delta t}+{\widehat{\mathrm{e}}}_0-{\widehat{\mathrm{v}}}_{\mathrm{mov}}=0\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill & {\widehat{\mathrm{i}}}_{L_0}+\frac{{\widehat{\mathrm{v}}}_{\mathrm{mov}}-{\widehat{\mathrm{e}}}_0}{R_0}-C\frac{e_0^{n+1}-e_0^n}{\Delta t}-{\widehat{\mathrm{i}}}_{A0}({\widehat{\mathrm{e}}}_0)-{\widehat{\mathrm{i}}}_{L_1}-\frac{L_1}{R_1}\frac{{\mathrm{i}}_{L_1}^{n+1}-{\mathrm{i}}_{L_1}^n}{\Delta t}=0\hfill \end{array}$$

(14b)

$$\begin{array}{cc}\hfill & {\widehat{\mathrm{i}}}_{L_1}+\frac{L_1}{R_1}\frac{{\mathrm{i}}_{L_1}^{n+1}-{\mathrm{i}}_{L_1}^n}{\Delta t}-{\widehat{\mathrm{i}}}_{A1}\left({\widehat{\mathrm{e}}}_0-L_1\frac{{\mathrm{i}}_{L_1}^{n+1}-{\mathrm{i}}_{L_1}^n}{\Delta t}\right)=0,\hfill \end{array}$$

(14c)

where we introduce the shorthand notation $\widehat{\mathrm{v}}=\frac{{\mathrm{v}}^{n+1}+{\mathrm{v}}^n}{2}$ for averaged quantities. The set of Equation (14) holds for a single MOV installed on a single phase conductor; in case of multiple MOVs, each MOV requires a single set. In Equation (14), V_mov is the voltage at the input terminals of the MOV, while the current i_mov entering into the MOV can be computed as ${\widehat{\mathrm{i}}}_{L_0}+\frac{L_0}{R_0}\frac{{\mathrm{i}}_{L_0}^{n+1}-{\mathrm{i}}_{L_0}^n}{\Delta t}.$

We must pay careful attention to the coupling of the IEEE MOV model with the MTL. Using the substitution theorem, the MOVs with known currents can be included into the MTL network through an ideal current source j placed at the tower node where MOVs are installed. Thus, the set of Equation (14) must be solved in conjunction with the MTL voltage equation at the kth node where MOVs are installed that reads

$$\frac{\mathrm{d}\mathbf{i}}{\mathrm{d}x}+{\mathbf{C}}^{\prime }\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}+\mathbf{j}\delta \left(x\right)=0.$$

(15)

As illustrated in Figure 5, where the case of an MTL with three phase conductors and two shield wires is reported, the ideal current source ${\mathbf{j}}_k$ at kth node can be expressed as

$$\mathbf{j}=\left[\begin{array}{c}{\mathrm{j}}_1\\ {\mathrm{j}}_2\\ {\mathrm{j}}_3\\ {\mathrm{j}}_4\\ {\mathrm{j}}_5\end{array}\right]=\left[\begin{array}{c}{\mathrm{i}}_{\mathrm{mov},A}\\ {\mathrm{i}}_{\mathrm{mov},B}\\ {\mathrm{i}}_{\mathrm{mov},C}\\ G_{\mathrm{s}1}\left({\mathrm{v}}_{k,4}-{\mathrm{v}}_0\right)\\ G_{\mathrm{s}2}\left({\mathrm{v}}_{k,5}-{\mathrm{v}}_0\right)\end{array}\right],$$

(16)

where

$${\mathrm{v}}_0=R_{\mathrm{g}}\sum _{m=1}^5{\mathrm{j}}_m$$

(17)

and $R_{\mathrm{g}}=G_{\mathrm{g}}^{-1}$ is the grounding resistance of the tower. Using a central difference combined with an average in time, we obtain the discrete version of Equation (15)

$$\left(\frac{{\mathbf{i}}_{k+\frac{1}{2}}^{n+1}-{\mathbf{i}}_{k-\frac{1}{2}}^{n+1}}{2}+\frac{{\mathbf{i}}_{k+\frac{1}{2}}^n-{\mathbf{i}}_{k-\frac{1}{2}}^n}{2}\right)+{\mathbf{C}}^{\prime }\Delta x\left(\frac{{\mathbf{v}}_k^{n+1}-{\mathbf{v}}_k^n}{\Delta t}\right)+\widehat{\mathbf{j}}=0.$$

(18)

The non linear system must be solved with three unknowns for each MOV (namely, $i_{{\mathrm{L}}_0}$, $i_{{\mathrm{L}}_1}$, and $e_0$) and ${\mathrm{v}}_0$.

Generally speaking, the winning approach for solving nonlinear circuits is to not solve their nonlinear equations directly but, instead, to linearize the circuits, build their corresponding linear equations, solve them, and repeat until convergence. In this case, we use a fourth-order extrapolation method to obtain an initial guess of the unknown variables with reasonable accuracy. Then, the nonlinear system in Equation (14) is solved through the NEQNF IMSL subroutine that uses a modified Powell hybrid algorithm and a finite-difference approximation to the Jacobian.

The exact procedure would require the solution of a matrix system of $N_w\times \left(2N_{\Delta \mathrm{x}}+1\right)$ linear equations along with $3\times N_{\mathrm{p}}\times \left(N_{\mathrm{span}}+1\right)$ non-linear equations, in the worst case, related to three surge arresters installed at each of the towers. However, generally, the surge arresters are installed only at a few (one or two, and a maximum of three) locations along the line. The inclusion of nonlinear elements makes the network equations to be a system of nonlinear equations. Solving such systems is not trivial and, in fact, is much harder than solving systems of linear equations. In addition, the tri-diagonal nature of the $\mathbf{A}$ matrix is lost along with the key advantages related to its solution.

We propose a multi-step strategy for simplifying the scheme and thereby reducing its computational cost. It consists of two steps. First, we solve Equations (12) and (16) for the unknowns ${\mathbf{v}}^{n+1}$ and ${\mathbf{i}}^{n+1}$ through a linear tri-diagonal system as in Equation (6). The matrix $\mathbf{A}$ is left unchanged regardless the possible presence of surge arresters; this is done by replacing the current generators ${\widehat{\mathbf{i}}}_{\mathrm{mov}}$ in Equation (16) with their approximation ${\mathbf{i}}_{\mathrm{mov}}^n$ at the previous time-step $t=n\Delta t$. Then, we use the current values ${\mathbf{i}}_{k\pm \frac{1}{2}}^{n+1}$ at the $N_{\mathrm{span}}+1$ tower nodes to solve the non-linear Equations (14) and (16) and compute the MOVs variables, update the ${\mathbf{v}}_k^{n+1}$ voltages at the tower nodes, and compute the current generators ${\mathbf{i}}_{\mathrm{mov}}^{n+1}$ to be used at the successive time-step. The system proved to be unconditionally stable and accurate, as shown in the next section.

3. Results

To validate and test the proposed CN-FDTD approach, we simulated a lightning first strike impinging on the shield wire of an MTL with nominal voltage of 132 kV. We analyzed the typical geometry depicted in Figure 1: it is composed of four wires, three phase conductor, and one shield wire. The radius of phase conductors is 15.75 mm, while that of the shield wire is 5.75 mm. The pole-to-pole span length ℓ is equal to 300 m and the overall length L of the line is assumed to be 9.3 km, namely 31 spans. The line is closed at both ends on its resistive dense characteristic impedance matrix [40]. The insulator gap length is 1.35 m.

The SW can be grounded at poles, as depicted in Figure 6. The maximum value of the tower footing resistance is typically below 20 $\Omega$: in our simulations, we assumed the towers grounded with a grounding resistance $R_{\mathrm{g}}=1/G_{\mathrm{g}}$ equal to $5\Omega$ or $10\Omega$. The MOVs are ABB PEXLIM P-G 120 kV.

We assumed that the strike hits the shield wire at distance $d=4.650$ $\mathrm{km}$ from the MTL endpoints, in the middle of the central span. The lightning current is a typical Heidler impulse with parameters ${\tau }_1=1.8$ $\mathsf{\mu }\mathrm{s}$, ${\tau }_2=95$ $\mathsf{\mu }\mathrm{s}$, $n=2$, $I_0=28$ $\mathrm{kA}$ [55,56]. The time waveform of the lightning current is shown in Figure 7. We did not take into account the presence of the AC voltages.

We investigated three different scenarios:

Scenario A:

The shield wire is grounded only on Tower #1 where a single MOV is installed on Phase Conductor A.

Scenario B:

The shield wire is grounded only on the first Tower #1 where two MOVs are installed on Phase Conductors A and B.

Scenario C:

The shield wire is grounded only on the first Tower #1 where three MOVs are installed on Phase Conductors A–C.

This simplified assumption avoided multiple reflections among towers due to the periodic grounding and allowed us to better focus on the role played by the combination of protection practices, i.e., SW grounding and MOVs.

Back-flashover occurs when the voltage across the insulator string—i.e., the difference between the voltage of the tower assumed equal to that of the SW ${\mathrm{v}}_4$ and the voltage of the phase conductor ${\mathrm{v}}_i$ with $i=1,2,3$—is higher than the electrical impulse withstand strength (or IL) of the insulator. In insulation coordination studies, the IL is typically assumed equal to the line Critical Flashover (CFO) multiplied by a factor 1.5, to consider the turn-up in the insulation volt–time curve for short front-time surges (standard waves are not longer than 50 $\mathsf{\mu }\mathrm{s}$).

The withstand strength of insulator strings is higher under surge voltages and pulses of short duration, while it reduces when the transient surges have longer duration. Without making use of probabilistic approaches and statistical procedures, a simplified expression of the withstand voltage capability for an insulator string is $V=\alpha +\beta /t^{0.75}$, with $\alpha =400L$ and $\beta =710L$, where L is the length of the insulator string and t is the elapsed time after lightning strike. The volt–time model for the insulator of the considered 132 kV-rated voltage line is reported in Figure 8. The volt–time model agrees well with the standardized test voltages suggested by IEC for testing equipment; for equipment with highest voltage equal to 145 kV, the standard lightning impulse withstand voltage (or IL) is equal to 650 kV and the reduced insulation level is 550 kV. In the adopted model, at 50 $\mathsf{\mu }\mathrm{s}$, the withstand voltage is 590 kV.

Figure 9 shows the voltage differences drops across insulators on Tower #1 when the line is protected only by the SW that is not grounded. In such a case, no MOV is installed. It is well known that in this case (that has little practical interest) a direct strike on the SW results in an insulation failure.

Figure 10 show the same voltage drops across insulators on Tower #1, when the line is protected only by the SW that is now grounded only at the first tower, with $R_{\mathrm{g}}=5$ $\Omega$ or $R_{\mathrm{g}}=10$ $\Omega$, respectively. MOVs are not installed yet. We observe that the transient overvoltages are deeply reduced; they are comparable in magnitude with the induced surge on lower Phase C and slightly higher than that on the highest Phase A. Obviously, the overvoltages increase linearly with increasing ground resistance $R_{\mathrm{g}}$. The overvoltages are not high enough to cause insulation failures, according to the proposed model, unless the ground resistance is increased above a certain value. However, they are high enough to cause the switching of MOVs.

Of special interest are the results in Figure 11, Figure 12 and Figure 13 that show the voltage drops across insulators on Tower #1 in the aforementioned Scenarios A–C, respectively. When the ground resistance $R_{\mathrm{g}}$ is low enough, e.g., 5 $\Omega$, the surge arresters switch for a limited time because the periodical grounding of the SWs is effective to mitigate the surge effects on the insulator strings. The overvoltages are limited below 100 kV by the insertion of MOVs. For higher values of the ground resistance, the grounding of the SWs reveals to be insufficient and the switching of the MOVs is clearly visible in the results. Interestingly, we do not see much value in increasing the number of MOVs from one to two; thereby, from a practical standpoint, it is preferable to use a single MOV per each tower or three MOVs, when a more reliable and effective protection is needed.

The CN-FDTD method proved to be robust and accurate. The results in Figure 11, Figure 12 and Figure 13 are reported for Courant numbers $f_{\mathrm{c}}$ equal to 1, 5, and 10. With the growth of the Courant number, we observe an accentuated dispersion of the scheme, clearly visible in the increasing time delay of the transient voltage. However, we do not see any particular difference in the absolute value for the transient overvoltages. The scheme is suitable for solving a large number of simulations (e.g., by means of Monte Carlo approaches) when performing a worst-case analysis in the design and coordination of line insulation.

In Figure 14, we consider a more practical case: the three towers are equipped with a single MOV on alternating phase conductor, i.e., we have a MOV on Phase A on Tower #1, a MOV on Phase B on Tower #2, and a MOV on Phase C on Tower #3. The grounding resistance is assumed to be 10 $\Omega$. Figure 14a shows the overvoltages on the insulators of the first tower. The overvoltages are more or less the same as in Figure 11b; the reflections due to the grounding of the SW on Towers #2 and #3 have little effects, because the grounding at the first tower mainly determines the protection. This is confirmed observing the overvoltages on the subsequent Towers #2 and #3 shown in Figure 14b,c: they are reduced to such a low value by the combination of SW grounding and MOV on the first tower that the surge arresters do not play any role.

In the same figures, we show the comparison between the CN-FDTD results and the predictions obtained through a classical leap-frog explicit FDTD scheme with $f_{\mathrm{c}}=1$. The agreement is excellent, validating the proposed implicit unconditionally stable scheme. The superiority of the CN-FDTD in terms of computational times has been deeply investigated [57]. A direct comparison of the explicit FDTD method with $f_{\mathrm{c}}=0.5$ and the Crank–Nicolson method with $f_{\mathrm{c}}=5$ indicates that less than half the number of computations produced a result with less than one-twentieth of the maximum error.

Finally, we consider a common situation in mountainous areas, where it is difficult to ground the SW at each tower. In this case, SW are grounded every three/four towers (i.e., three spans) with $R_{\mathrm{g}}=10$ $\Omega$ and one or three MOVs are installed on the same towers where SWs are connected to the Earth. As shown in Figure 15, we assume that the SWs are grounded at Towers #1 and #4 and the lightning strikes in the midspan between Towers #2 and #3. Figure 16 shows the overvoltages across insulators on Tower #3 (tower without MOVs and grounding system) and Tower #4 (tower with MOVs and grounding system) when one MOV is installed on Phase A or three MOVs are installed on all phase conductors on Towers #1 and #4. Multiple reflections of the voltage waves at the grounded towers are clearly visible. Comparison of the volt–time characteristic in Figure 8 with the peak values of the surges on Tower #1 clearly shows that the insulation fails due to a back-flashover. The installation of one or three MOVs on Tower #2 is only partly beneficial on Tower #1; they are effective only to improve the insulation level on Tower #2.

4. Conclusions

Simulation of metal oxide surge arrester is a challenge for the proper computation of transients in overhead power lines. Several factors make metal oxide varistor modeling a difficult task: their behavior is non-linear and strongly frequency dependent; many topological variations on materials and construction are possible; there are many physical and parasitic attributes whose behavior may need to be correctly represented; and several MO surge arresters present hysteresis-type loop, because the discharge voltage reaches the peak before the discharge current reaches its own peak. Various popular models have been proposed during the last years to reproduce the performance of an MO surge arrester in high-frequency transients. Additionally, the inclusion of a surge arrester model into the multiconductor transmission line is not a straightforward task and several additional aspects should be considered for an accurate simulation.

This paper proposes an efficient simulation of the MOV model recommended by IEEE, based on a smart reformulation of system nonlinear equations that allows to achieve a faster solution convergence, while maintaining stability and high accuracy. The model allows simulating the role played by the tower grounding resistance, and, in general, impedance on the MOV performance, and, consequently, conducting a coordinated study on the effects of ground potential rise and MOVs on insulation level to assess the lightning performance of the line.

The model was inserted into a new unconditionally stable Crank–Nicolson finite difference time domain scheme that was developed for the simulation of transient propagation over overhead lines; it ensures the exceptional robustness and accuracy. The approach is not without its pitfalls, since it has matrix computational complexities to be dealt with. Nevertheless, the symmetric diagonally dominant block tri-diagonal form of the solving matrix leads to a fast solution enhanced by the use of modern GPU parallel solvers.

As noted above, the implicit scheme not only allows the reduction of the total running time through the choice of larger Courant factors, but also preserves its stability despite high-order current and voltage harmonics may rise due to non-linearities.

With reference to the configuration shown in Figure 6, the running time on a Intel workstation of the Yee-based explicit scheme with $f_{\mathrm{c}}=0.8$ was about $13.5$ s, while 8 and 5 s were the times required for the simulations carried out by the Crank–Nicolson implicit scheme with $f_{\mathrm{c}}=5$ and $f_{\mathrm{c}}=10$, respectively. Thereby, a reduction of the running time of up to $62\%$ is achieved thanks to the relaxed time discretization allowed by the implicit method.

As a consequence, the numerical implicit scheme proposed in this paper allows for assessing the beneficial effects of MOVs to the HV insulation coordination and may further prove its convenience when a Monte Carlo procedure is applied to compute probability distribution functions for key outputs.