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

: Surge arresters may represent an efﬁcient choice for limiting lightning surge effects, signiﬁcantly reducing the outage rate of power lines. The present work ﬁrstly presents an efﬁcient numerical approach suitable for insulation coordination studies based on an implicit Crank–Nicolson ﬁnite 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 efﬁcient 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 ﬁnite 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).


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 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.

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): where v and i are the voltages and currents at coordinate x along the line and at time instant t, respectively; L ext is the external line inductance per unit length (p.u.l.); Z ζ 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 C is the p.u.l. geometric capacitance. The matrix C can be replaced with a non-linear time-dependent dynamic capacitance C 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 G, since transversal losses are usually associated with the corona effect and accounted for through the dynamic capacitance C 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 w wires (see Figure 1), i.e., N p phase conductors and N s overhead shield wires (OHSWs), with N w = N p + N s . The length of the line is L = N span , where is the pole-to-pole span and N span is the number of sections. The OHSWs are grounded at some poles by a constant equivalent ground admittance G 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 e,i of the wire computed as h e,i = h i − 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 µ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.

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 ∆x, where we seek voltages v k at N ∆x + 1 integer mesh points x k (the mesh boundaries are voltage nodes) and currents i k+ 1 2 at N ∆x half-integer mesh points x k+ 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 ∆t into two sub-steps ∆t/2. In the first sub-step the explicit Forward Euler method is used to deduce the unknowns at time n + 1 2 ∆t: In the second sub-step, we use the implicit Backward Euler method to compute the unknowns at time (n + 1) ∆t: 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 i n+1 The last two equations (which are valid for lossless MTLs) make evident that the CN scheme is based on a central difference on time n + 1 2 ∆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: We solve Equation (5) alternating a voltage node with a current node in the unknown column vector x n+1 and we obtain a linear system in the form: In Equation (6), A is a block tri-diagonal matrix (the size of each block is N w × N w ) whose band is 2N w − 1: If the p.u.l. parameters are time-independent (i.e., we neglect corona effect that results in the time-dependent dynamic capacitance matrix C dyn (t)), the matrix 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 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 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 n + 1 2 ∆t as If we use a trapezoidal rule to approximate the second integral over the last time-step and we denote the first history integral as H(x, t), we end up with 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: 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:

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 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.
R 0 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 bp ) and stray capacitance (C 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: 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 µ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 µ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 L 0 and i 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: 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): where we introduce the shorthand notationv = v n+1 +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î 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) 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 j k at kth node can be expressed as where and R g = G −1 g 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)   i n+1 The non linear system must be solved with three unknowns for each MOV (namely, i L 0 , i L 1 , and e 0 ) and 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 × (2N ∆x + 1) linear equations along with 3 × N p × N span + 1 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 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 v n+1 and i n+1 through a linear tri-diagonal system as in Equation (6). The matrix A is left unchanged regardless the possible presence of surge arresters; this is done by replacing the current generatorŝ i mov in Equation (16) with their approximation i n mov at the previous time-step t = n∆t. Then, we use the current values i n+1 k± 1 2 at the N span + 1 tower nodes to solve the non-linear Equations (14) and (16) and compute the MOVs variables, update the v n+1 k voltages at the tower nodes, and compute the current generators i n+1 mov to be used at the successive time-step. The system proved to be unconditionally stable and accurate, as shown in the next section.

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 Ω: in our simulations, we assumed the towers grounded with a grounding resistance R g = 1/G g equal to 5 Ω or 10 Ω. The MOVs are ABB PEXLIM P-G 120 kV.
We assumed that the strike hits the shield wire at distance d = 4.650 km from the MTL endpoints, in the middle of the central span. The lightning current is a typical Heidler impulse with parameters τ 1 = 1.8 µs, τ 2 = 95 µs, n = 2, I 0 = 28 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.  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 v 4 and the voltage of the phase conductor 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 µ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 = α + β/t 0.75 , with α = 400L and β = 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 µ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 g = 5 Ω or R g = 10 Ω, 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 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 Figures 11-13 that show the voltage drops across insulators on Tower #1 in the aforementioned Scenarios A-C, respectively. When the ground resistance R g is low enough, e.g., 5 Ω, 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.     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 Ω. 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 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 c = 0.5 and the Crank-Nicolson method with f 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 g = 10 Ω 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.

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 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 c = 5 and f 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.