Monte Carlo method for estimating backflashover rates on high voltage transmission lines

Electric Power Systems Research 119 (2015) 247–257

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

Monte Carlo method for estimating backflashover rates on high voltage transmission lines Petar Sarajcev ∗ University of Split, FESB, Department of Power Engineering, R. Boskovica 32, HR-21000 Split, Croatia

a r t i c l e

i n f o

Article history: Received 12 May 2014 Received in revised form 6 October 2014 Accepted 9 October 2014 Keywords: Backflashover Transmission line Monte Carlo method EMTP-ATP Bivariate log-normal distribution Insulator string flashover

a b s t r a c t This paper presents a novel Monte-Carlo based model for the analysis of backflashover rate (BFOR) on high voltage transmission lines. The proposed model aims to take into the account following aspects of the BFOR phenomenon: transmission line (TL) route keraunic level(s), statistical depiction of lightningcurrent parameters (including statistical correlation), electrogeometric model of lightning attachment, frequency-dependence of TL parameters and electromagnetic coupling effects, tower geometry and surge impedance, tower grounding impulse impedance (with soil ionization), lightning-surge reflections from adjacent towers, non-linearity of the insulator strings flashover characteristic, distribution of lightning strokes along the TL span and power frequency voltage. In the analysis of the BFOR, special attention is given to the influences emanating from the insulator strings flashover characteristic and lightning statistics. The model could be applied to the transmission line as a whole or some of its portions, e.g. first several towers emanating from the substation or several towers crossing a mountain ridge. © 2014 Elsevier B.V. All rights reserved.

1. Introduction High voltage (HV) transmission lines are exposed to lightning strikes, where only direct lightning strikes (to shield wire(s), phase conductors and tower tops) are of engineering concern; nearby lightning strikes have no influence on the HV transmission line (TL) performance. Direct lightning strikes to phase conductors, where shield wire(s) is(are) present on the tower, are accompanied by lightning currents with constrained amplitudes, due to the shielding effect of wire(s); nonetheless, these strikes can provoke a flashover of the TL insulation (i.e. insulator string flashover). The rate at which this is to be expected, per 100 km-years of transmission line, is termed the shielding failure flashover rate (SFFOR). Direct lightning strikes to shield wire(s) and tower tops can also provoke a flashover of the TL insulation, where the strikes to the tower tops are more significant in producing insulator flashovers (and statistically speaking more probable) then the strikes to midspans. The rate at which this is to be expected, per 100 km-years of transmission line, is termed the backflashover rate (BFOR). The flashovers away from insulator strings, both in SFFOR and BFOR analysis, are regarded as being far less probable than the insulator string flashovers; hence, only flashovers at the TL insulator strings are considered possible in numerical treatment of the

∗ Tel.: +385 21 305806; fax: +385 21 305776. E-mail address: [email protected] http://dx.doi.org/10.1016/j.epsr.2014.10.010 0378-7796/© 2014 Elsevier B.V. All rights reserved.

phenomenon. It is the intention of this paper to analyse solely the BFOR of HV transmission lines. The backflashover occurrence rate is important, along with the mentioned SFFOR, for estimating the outage times/rates of transmission lines due to lightning. It is important in designing the HV substations (or switchyards) overvoltage protection, in terms of the incoming overvoltage emanating from the backflashovers on neighbouring TL towers incident to the station. Furthermore, it is of importance in the decision making process regarding the shielding of TLs using surge arresters (TLA applications on specific parts of the TL route). It has been extensively studied by many researchers, using analytical and numerical methods, and the volume of published material on the subject is overwhelming. The analytical methods are extensively described by the IEEE WGs [1,2] and CIGRE WGs [3], with additional details provided in numerous references cited therein. A comparison between these recommendations in given in Ref. [4]. Further extensive exposition of analytical methods is provided in Refs. [5, Ch. 10]. Nowadays, it is far-more common to treat the backflashovers on TLs in terms of the numerical simulations, carried-out by means of the Electromagnetic Transients Programs (EMTP), e.g. [6–8]. With the numerical approach to the transient analysis of TL lightning surges, detailed (and often quite sophisticated) models of the TL components are needed, some of which exhibit non-linear behaviour, frequencydependence, etc. The IEEE WGs and CIGRE WGs offer extensive guidelines when it comes to representing transmission line elements (and other network elements) for numerically simulating

248

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

fast-front (i.e. high-frequency) transients [9,10]. Furthermore, special recommendations exist for particular network components. Interested reader is at this point advised to consult the extensive treatment of modelling guidelines for TL lightning-surge numerical simulations provided in [11, Ch. 2] and references cited therein. Further important simulation details, concerning the backflashover analysis on HV transmission lines, can be found in Refs. [12–17]. Direct lightning strike to the transmission line tower (or its near-vicinity) initiates a travelling wave process in the system of lightning channel, shield wires, TL tower, and tower footing impedance. Accompanying current and voltage transient states, at the various points of the TL tower, are established through the complex propagation pattern of current and voltage travelling waves, including reflection and transmission of those waves on various points of travelling paths wave impedance discontinuities. This process results with a transient voltage being applied on the TL insulator strings with the possibility of their flashover, i.e., when the transient potential of the tower arm exceeds the critical flashover voltage of the insulator strings (itself being a non-linear function of the applied voltage), biased to some extent by the power-frequency phase voltage, it causes a backflashover. The numerical simulation of the TL backflashover proposed in this paper aims to take into the account following aspects of the phenomenon: transmission line route keraunic level(s), statistical depiction of lightning-current parameters (including statistical correlation between the parameters), electrogeometric model (EGM) of the lightning attachment process, frequency-dependence of TL parameters and electromagnetic coupling effects, tower geometry and surge impedance, tower grounding impulse impedance (with soil ionization if present), lightning-surge reflections from adjacent towers, non-linear behaviour of the insulator strings flashover characteristic, TL span length, statistical distribution of lightning strokes along the TL span and power frequency voltage. Furthermore, the proposed approach aims at coupling the Monte Carlo method with the EMTP simulation, in order to establish the statistical probability of backflashovers on HV transmission lines. The Monte Carlo method has been applied to the problem of calculating BFOR, for example, in Refs. [12,18–22]. However, the here proposed approach is unique in the way it implements the Monte Carlo method and in the way in which it treats the statistical parameters of lightning currents incident to TLs (accounting for the electrogeometric model of lightning attachment along with statistical dependence between lightning-current parameters). On top of that it implements a state-of-the-art TL model for the EMTP backflashover simulation. Monte Carlo procedure applied, for example, in Ref. [22] generates lightning data from the downward (negative) lightning statistics, chooses at random lightning starting points up to some distance from the TL, and, using EGM (from Eriksson), determines whether this lightning will strike the TL or the nearby earth. If there is a strike to the TL, the EMTP simulation is carried-out. Here proposed method, on the other hand, first computes (numerically) the probability density function of the statistical distribution of lightning currents which are, by means of applying the EGM, incident to TL and then from it generates lightning data (accounting for the statistical correlation between parameters) for the EMTP simulations. This reduces the number of samples used for the analysis. The proposed method, through simulations, provides insight into the BFOR behaviour—in the statistical sense—which is due to many influential factors, some of which can be changed between simulation runs. Furthermore, considering the well-established influence of the tower grounding impulse impedance and insulator strings flashover characteristic on the BFOR, these aspects will be numerically investigated within the sensitivity analysis provided in the paper. The paper is organised in the following manner. In Section 2, a brief outline of the TL model for the BFOR analysis is provided,

which is suitable for the implementation in the EMTP-ATP software package. Section 3 provides necessary statistical treatment of the lightning current parameters. In Section 4 is provided an estimation of the number of direct lightning strikes to transmission line, by means of implementing the electrogeometric model of lightning attachment in combination with the statistical distribution of lightning strokes along the TL span length. This section additionally presents the probability density functions of lightning current amplitudes incident to transmission lines, in accordance with the EGM attachment process. Section 5 provides the information on the implementation of the Monte Carlo method and its coupling with the EMTP simulation running in the batch mode. A test case of the HV transmission line, along with the sensitivity analysis, is provided in Section 6, which is followed with the conclusion in Section 7. 2. Transmission line modelling for backflashover analysis The EMTP model of the HV transmission line for lightning surge transient simulation in general, and backflashover analysis in particular, has been thoroughly studied and widely published, see Refs. [8–10,23,24]. A brief outline of the EMTP-ATP model, as employed for the purpose of this paper, will be presented in this Section. The model consists of several components: (i) TL phase conductors and shield wire(s), including spans, line terminations and power frequency voltage, (ii) TL tower, (iii) tower grounding impedance, (iv) insulator string (i.e. archorn) flashover characteristic, (v) lightning current and lightning-channel impedance. A screenshot of the main part of the TL model, constructed within the ATPDraw preprocessor to the EMTP-ATP software package, is presented in Fig. 1. 2.1. Phase conductors, shield wire(s), spans, line termination, power frequency voltage High voltage transmission line phase conductors and shield wire(s) are modelled as distributed-parameters, untransposed, frequency-dependent, multiphase transmission line, by means of employing the so-called LCC component of the EMTP-ATP which utilises the JMarti TL model [6,7,25]. Phase conductors and shield wire(s) positions on the tower (from the most-representative tower within the TL route) are used, along with their maximum allowed sags, cross-sectional dimensions, DC resistances, ground resistivity of the ground return path, etc. The electromagnetic coupling effects of TL phase conductors and shield wire(s) are accounted for by this model. However, the effects of corona on the propagation of lightning surge on the conductors are neglected in this model. Five spans of the transmission line, at each side of the tower being struck by lightning, are modelled in this way, using eight decades with ten points per decade in the JMarti frequency-dependent model, with the modal transformation matrix computed at the (dominant) frequency of 400 kHz [6,7]. After the fifth span, at each side of the struck tower, additional 30 km of the transmission line is modelled in the same way, and the line model is then terminated by an ideal, grounded, power-frequency, three-phase voltage source. 2.2. Transmission line towers The steel-lattice towers of HV transmission lines are usually represented as a single conductor, distributed-parameter (frequency-independent) lines. The towers of EHV and UHV transmission lines, on the other hand, are represented using more complicated so-called multi-story structures, which are not considered adequate for HV transmission lines, e.g. [14]. The single value of the tower surge impedance is computed from the analytical expressions, based on the theoretical background provided in [26], which depend on the tower configuration and can be found

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

249

lightning strike, where reflections from the tower base will arrive much sooner at the tower top then reflections from adjacent towers. Hence, the influence of the (apparent) TL tower footing (i.e. grounding) surge impedance on the tower top transient voltage is determined by its response time and current dependence. In cases where the tower grounding (system) cannot be regarded as being concentrated (e.g. there are counterpoises installed covering distances greater than cca. 30 m), frequency-dependent tower grounding impedance model is needed. Its implementation can be somewhat complicated [16]. On the other hand, within some 30 m of the tower base (i.e. concentrated tower grounding system), the tower grounding impedance exhibits only current dependence, which can be modelled in accordance with guidelines provided in [9] and implemented in EMTP-ATP software package by means of the MODELS language [25]. The following equation is utilised for that purpose [9]: Rtower =

R0



1 + It /Ig

(1)

where R0 is the tower grounding resistance at low frequency and low current magnitudes (), It is the lightning current through the tower footing impedance (kA), and Ig is the lightning current level which determines the soil ionization inception process (kA). This current is determined from the following expression [9]: Ig =

E0 2R02

(2)

where E0 is the soil ionization electric field gradient, provided in (kV/m) and  is the apparent soil resistivity in (m). It is evident that the inception of the soil ionization depends on the value of the so-called soil ionization electric field gradient and the low frequency value of the tower grounding resistance (which in-turn depends on the soil resistivity and type of the grounding system employed). The value of the soil ionization electric field gradient usually ranges between 300 (kV/m) and 1000 (kV/m), with value of 400 (kV/m) being selected for the purpose of this paper. Grounding impedances of the adjacent towers (five of them on each side of the struck tower) are modelled with a simple resistance, the value of which is equal to the tower grounding resistance at low frequency and low current magnitudes. Namely, the influence of the models of these impedances on the struck tower’s top overvoltage formation is negligible. 2.4. Insulator strings flashover characteristic Fig. 1. Screenshot of the main part of the TL model constructed within the ATPDraw pre-processor to the EMTP-ATP software package.

in, e.g., Ref. [9]. The velocity of the surge propagation along the steel-lattice tower is assumed to be equal to the speed of light in free space (although some authors assume somewhat lower value). Additionally, tower arms could be modelled as transmission line stubs (distributed parameter, single conductor) with constant value of surge impedance; some authors use inductances instead. However, the influence of the tower arms could be neglected altogether. Each of the towers, eleven of them in total, (struck tower plus five towers at each side) is modelled in this way. Tower tops are connected to shield wires. 2.3. Tower grounding impedance The model of the lightning-struck TL tower grounding impedance is important factor influencing the subsequent formation of the overvoltage on its tower top (and its arms), due to subsequent reflections of travelling waves formed by the said

The insulator strings (i.e. archorn) flashover characteristic is a non-linear function of the applied impulse voltage and it exhibits complicated behaviour in nature, which is rather difficult to fully reproduce. Interested reader is advised to consult Refs. [13,27–29] for more information. Hence it is usually modelled in the EMTPATP software package by means of the voltage-controlled switches. The flashover characteristic itself is programmed using the MODELS language [25]. Sometimes, a parallel capacitor is added to the controlled switch in order to simulate the coupling effects between conductors and tower structure [9]. The model of the insulator strings flashover has an important effect on the backflashover analysis (i.e. its probability of occurrence), which will be thoroughly investigated in this paper. The following three different models of the insulator strings flashover characteristics will be implemented in the backflashover analysis provided within the context of this paper: (1) voltage-time characteristic, (2) equal-area criterion, and (3) leader progression method. Each of these models are implemented within the MODELS language and used to control a TACS switch [25]; closure of this switch represents a flashover of the insulator strings.

250

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

The voltage-time characteristic of the insulator strings can be depicted by the following analytical expression [9]:

v(t) = 400 · s +

700 · s in (kV) t 0.75

(3)

where s is the length of the insulator strings (i.e. length of the insulator archorns gap) in meters, and t is the elapsed time after the lightning stroke in (␮s). The equal-area criterion in determining the insulator flashover characteristic, also termed the integration method, utilises following equation [9]:



tc n (v(t) − V0 ) · dt

D=

(4)

However, due to the fact that lightning is stochastic in nature, above mentioned parameters can only be described in statistical terms. In HV transmission line studies involving lightning strikes, including the backflashover analysis, only negative downward lightning strikes are of engineering interest. Hence, only parameters of this lightning current type will be presented hereafter. These parameters are provided, among others, in the following publications [31–35]. It has been established that probability distributions of the lightning-current parameters (amplitudes, front and tail durations) individually follow a log-normal distribution, in which case the probability density function (PDF) of the random variable can be given by the following expression [33]:

t0

where t0 is the time after which the voltage v(t) is higher that the required minimum voltage V0 , also known as the reference voltage; tc is the time-to-breakdown, and D is the disruptive effect constant. If n = 1 than the method is also known as the equal-area law. The difficulty with the application of this method is in choosing the right values of D, V0 for the insulation strings being investigated. The value of the integral in (4) is monitored during the EMTP-ATP simulation, by means of the MODELS language, and if it attains or exceeds the value of D, the associated TACS switch is closed, signifying the occurrence of the insulator strings flashover. The most-sophisticated model of the insulator strings flashover behaviour is the so-called leader development or leader progression model. There are several variants of this model, e.g. see Ref. [30], but the particular one used in this paper is based on the solution of the following differential equation [9]: d = k · v(t) · dt



v(t) dg −  (t)



− E0

f (x) =

2 ] exp[−(ln x − ln )2 /2ln x √ 2xln x

where  represents the median value and  lnx represents the associated standard deviation of the ln x. Each of the three lightning-current parameters can be individually depicted by a lognormal distribution with appropriate parameters (median value and standard deviation). However, situation is complicated by the fact that there has been found a statistically significant correlation between the lightning-current amplitudes and front durations. This necessitates usage of the joint (i.e. bivariate), as well as conditional, probability distributions in their treatment. The joint (i.e. bivariate) probability density function, in case of the lightning current amplitude (I) and front duration (tf ), can be described by the following relation [33]: f (I, tf ) =

(5)

where dg is the insulator strings (i.e. archorn) length,  (t) is the leader length, v(t) is the actual (absolute value) voltage on the insulator strings, and k, E0 are constants which are found to be dependent on the type of the insulator. The differential equation is solved during the EMTP-ATP simulation, by means of the MODELS language, for the length of the leader at each simulation time-step. If this length attains or exceeds the gap length (i.e. the length of the insulator strings), the associated TACS switch is closed, signifying the occurrence of the insulator strings flashover. 2.5. Lightning current A direct lightning strike to the TL tower top is, in the EMTPATP software package, usually modelled as an ideal current source in parallel with the resistance, the value of which represents the lightning-channel surge impedance. This value ranges between 400 (m) and 1000 (m) and is hereafter assumed at 400 (m). The current source can be of different types: ramp-slope type, CIGRE type, double-exponential type, Heidler type, to name the few most-prominent ones [25]. Using a different type of the current source will result with somewhat different overvoltages being formed on the TL tower top and associated elements. The Heidler type of lightning-current source will be utilised for the purpose of this paper, due to the fact that it is generally recommended for lightning-associate transient simulations. This ideal current source is connected to the top of the TL tower being struck by lightning (i.e. middle tower of the TL model). 3. Statistical depiction of lightning current parameters In terms of the EMTP lightning-surge analysis, lightning current is depicted (along with a polarity), with the following three parameters: (1) amplitude, (2) front duration and (3) tail duration.

(6)

exp[−(f1 − f2 + f3 )/(2 · (1 − c2 ))] 2 · I · tf · ln I ln tf ·

with



f1 =

ln I − ln I ln I



f2 = 2c ·

 f3 =

1 − c2

(7)

2 (8)

ln I − ln I ln I

ln tf − ln tf ln tf





ln tf − ln tf ln tf

 (9)

2 (10)

where I ,  lnI represent median value and standard deviation of the lightning current amplitudes, tf , ln tf represent median value and standard deviation of the lightning current front durations, and c is the coefficient of correlation between the lightning current amplitudes and front durations. If the statistical variables are independently distributed, which is the case with the lightning current amplitude (I) and tail duration (th ), then associated c = 0 and Eq. (7) reduces to f(I, th ) = f(I) · f(th ), with f(I) and f(th ) obtained from (6) by introducing relevant median values and standard deviations. Further statistical treatment of lightning current parameters is presented in e.g. Ref. [33]. Following parameters for the statistical distributions of (negative downward) lightning current parameters will be utilised (hereafter termed the original set): I = 31.1 (kA),  lnI = 0.484; tf = 3.83 (␮s), ln tf = 0.55; c (I, tf ) = 0.47; th = 77.5 (␮s), ln th = 0.58; c (I, th ) = 0. These parameters are recommended in [33]. As an alternative, following parameters are provided for the log-normal distribution of lightning current amplitudes (hereafter termed the alternative set): I = 30.1 (kA),  lnI = 0.76 and for the front duration tf = 2.0 (␮s), ln tf = 0.494, with tail duration parameters and correlation coefficients inherited from the original set. This is due to the fact that there are differences between lightning-current parameters provided by different researchers, e.g., statistical parameters

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

of lightning strokes in Japan are quite different from those just provided.

The maximum shielding failure current can be determined from the following expression: Im =

4. Estimation of the number of direct lightning strikes to transmission line

r

gm

1/b (12)

A

with [5, Ch. 7,11, Ch. 6]:

The problem of estimation of the number of direct lightning strikes to transmission lines has been tackled by means of the electrogeometric model (EGM) of lightning attachment for quite some time, and the literature on the subject is exhaustive. The interested reader is at this point advised to consult Refs. [5, Ch. 7,36, Ch. 4], as well as, for example, [19,37–40], for the in-depth treatment of this subject. Only a brief outline of the methodology, as far as it is needed for the purpose of this paper, will be presented here. According to the EGM of lightning attachment to transmission lines, and considering only vertical lightning strikes, the number of lightning strikes to phase conductors and shield wire(s) depend on their exposure areas, which are determined in terms of the lightning striking distance and tower geometry. According to theory presented in Ref. [5, Ch. 7] following expression for estimating the number of direct lightning strikes to shield wire(s) can be obtained:

rgm =

(h + y)/2

NG = 2LNg ·



Im

∞

Dg (I)f (I) dI + 2LNg · 0



D g (I)f (I) dI + LNg Sg Im

(11)

where Ng = 0.04Td1.25 in (km−2 year−1 ) is the annual average ground flash density (Td is the long-term average annual number of thunderstorm days); f(I) is the probability density function of the lightning current amplitudes distribution; L is the transmission line length; Sg is the distance between shield wires (in case of a single shield wire Sg = 0); Im is the maximum shielding failure current; Dg (I) and D g (I) are exposure distances for the shield wire(s) as a function of lightning current amplitudes. The maximum shielding failure current and exposure distances for the shield wire(s) depend on the EGM that is being applied to the TL geometry (there are several possible EGM models to choose from). Strictly speaking, expression (11) is valid for TLs with horizontal conductor configurations but could be easily adopted for the TLs with vertical conductor configurations. If the TL route traverses through terrains with different keraunic levels, the route is then split into sections having different keraunic levels and Eq. (11) is solved for each of the sections separately. Lightning strikes are not uniformly distributed along the TL span length, which has been shown in Ref. [41]. If one considers only lightning strikes to TL tower tops (and their nearvicinity), which is often the case with the BFOR analysis, then the expression (11) needs to be corrected with the appropriate coefficient which takes into the account the actual statistical distribution of lightning strikes along the TL span length. For example, according to the analysis provided in Ref. [41] following expression is valid for HV transmission lines (above 230 kV): Pt = 5.486 × 10−7  2 − 7.452 × 10−4  +0.3587, where  is the length of the TL span in meters; for HV transmission lines below 230 kV the Pt value needs to be increased by a factor of 0.02. If one considers lightning strikes collected by the TL tower top and the shield wire(s) up to the 55% of the half-span from each side of the stricken tower, then, in accordance with [41], it yields that Pt = 0.6. The EGM brings into functional relationship the lightning current amplitude with its striking distance to phase conductors, shield wire(s) and to ground surface, and is of the form r = A · Ib (m), where A and b are model parameters and I is the lightning-current amplitude in (kA).

(13)

1 − (rc /rg ) · sin(tan−1 (a/h − y))

where rc and rg depict striking distances to phase conductors and ground surface, respectively; h is the height of the shield wire(s) on the tower; y is the height of the phase conductor(s) at the tower and a is the length of the tower arm(s) carrying the phase conductor(s), all supplied in meters. Exposure distances for the shield wire(s) can be, in accordance with the EGM model, determined from the following expressions [5, Ch. 7]:



Dg = rc · cos

a − sin−1 tan−1 h−y



rc2 − (rg − h)



Dg =

2



a2 + (h − y) 2rc

2



for rg ≥ h

(14)

(15)

for rg < h

rc



251

Introducing different EGMs in (11) and solving the associated integrals numerically yields the expected number of direct lightning strikes to transmission line, for any value of the average ground flash density, which can be further biased by the statistical distribution of lightning strikes along the TL spans (using the Pt factor introduced previously). Moreover, this same application of the EGM to the exact TL tower geometry yields the probability density function (PDF) of lightningcurrent amplitudes that are incident to the HV transmission line under scrutiny. This new PDF is determined from the following equation [5, Ch. 7]:

g(I) =

⎧ LN (2D + S ) ⎪ ⎨ g g g · f (I) for I ≤ Im NG

 ⎪ ⎩ LN g (2D g + Sg ) · f (I) for I ≥ Im

(16)

NG

in which f(I) stands for the PDF of the general (negative downward) lightning current amplitudes. It can be shown that this new PDF function g(I) depicts also a log-normal distribution. This can be proven through fitting the statistical distribution to the random data generated by (16), by means of applying the Maximum Likelihood Estimation of its parameters, e.g. [42, Ch. 2,43, Ch. 6]. Also, selection of the EGM model, which has been applied to the TL under scrutiny, influences to some extent the parameters of this new log-normal probability distribution function. Further details can be found in, e.g., Ref. [39]. Hence, lightning current amplitudes that are drawn from this new log-normal distribution g(I) are, by inference, incident to the TL under scrutiny, due to the fact that the associated PDF, depicted by (16), has been formed from the application of the EGM to the tower geometry of the concrete TL under investigation. 5. Monte Carlo method for the BFOR calculation The BFOR calculation procedure, utilising the Monte Carlo method, can be decomposed into three stages: pre-processing, numerical simulation and post-processing, with first and third stages implemented in a purposefully developed computer program. In the first stage, random variates representing lightningcurrent parameters are drawn from the univariate and bivariate log-normal statistical distributions. The second stage is implemented by means of the EMTP-ATP software package running

252

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

in batch mode, with interventions on its input and output files carried-out during the simulation, through the developed computer program. In the third stage, computation of the BFOR probability and BFOR per 100 km-years is performed, using the data provided by the numerical simulations and additional computations of the number of lightning strikes to transmission line. Drawing a random variate from the Log − N(x ,  lnx ) is accomplished by transforming variate from the standardised normal distribution, itself obtained with the Box–Muller method. However, when it comes to drawing a random variate from the modified PDF provided by (16), the Acceptance–Rejection method has been chosen due to its simplicity. In the case of generating correlated lightning-current parameters one needs to draw these random variates from the appropriate bivariate log-normal distribution, since they are statistically dependent. This is accomplished by means of transforming variates drawn from the standardised bivariate normal distribution [43, Ch. 4]. Namely, if the twodimensional statistical variable Y = [Y1 , Y2 ]T is drawn from the standardised bivariate normal distribution Y ∼ N(, ˙) where  = [0, 0]T is the mean vector and ˙ is the variance-covariance matrix



˙=

1

c

c

1



(17)

then the associated two-dimensional statistical variable from the bivariate log-normal distribution (denoting a lightning current amplitude and a front duration) can be determined as follows:

 X=

I · exp(ln I · Y1 ) tf · exp(ln tf · Y2 )

 (18)

where I , tf are median values of lightning current amplitudes and front durations;  lnI , ln tf are their standard deviations, while the variable c stands for the correlation coefficient between these statistical variables, in accordance with expressions (7)–(10). Preservation of the correlation coefficient between standardised bivariate normal and the appropriate bivariate log-normal distribution is tested-for by means of the Spearman’s correlation coefficient. In the case in which the statistical variables of the bivariate log-normal distribution are mutually independent (as is the case for the lightning-current amplitudes and tail durations) the drawing of random variates proceeds again by means of transforming random variates drawn from the standardised bivariate normal distribution, with the coefficient of correlation equal to zero in the associated variance-covariance matrix [43, Ch. 4]. A large number (N) of randomly generated lightning-current parameters are drawn from the appropriate statistical distributions in a manner described above, forming three random variates (amplitude, front duration and tail duration), with statistical correlation being present between amplitudes and front durations. Furthermore, in order to account for the power frequency phase voltages, phase angle of the three-phase ideal current source is randomly chosen for each simulation from the uniform distribution between 0◦ and 360◦ . For each set of lightning-current parameters (and phase voltages) an EMTP-ATP simulation of the backflashover is performed. For each simulation run of the EMTP-ATP a new set of lightning-current parameters (and phase voltage angles) need to be written into the associated. DAT file (initially created by the ATPDraw), then ATP file needs to be executed, and backflashover occurrence checked-for in the newly created. lis file [25]. This necessitates broadcasting position (open/closed) of the TACS switch associated with the insulator string flashover models, created by means of the MODELS language. If the backflashover occurs on any of the phases a BFOR counter for that phase is increased and

a cumulative sum is kept for each phase. This process is automated using the purposefully developed computer program. After all N simulations are performed, there is some k number of backflashovers for each of the phases. Hence, probability of backflashover, for some phase, is derived from the frequency of the backflashovers for that phase, i.e., p = k/N. Assuming Binomial distribution for this probability, a two-sided confidence interval for the probability can be determined from the following asymptotic relationship [42, Ch. 2]: pu po





=

1 N + 2q

k+

2q 2



∓ q

2q k(N − k) + N 4



(19)

where pu , po are lower and upper confidence limits, respectively, and q is the quantile of the standardised normal distribution of order q = (1 + )/2, with being the desired confidence coefficient. The BFOR for the TL under scrutiny can finally be estimated using the following relationship: BFOR = NG · [p ∓ (pu , po )]

(20)

with NG signifying the estimated number of direct lightning strikes into the transmission line, in accordance with (11), corrected for the influence of the statistical distribution of lightning strikes along the TL span length (by means of the coefficient Pt ). 6. Test case transmission line and sensitivity analysis Heretofore presented methodology for estimating the BFOR on HV transmission lines will be demonstrated on the typical single-circuit 110 kV transmission line with vertical conductor configuration and steel-lattice towers. Tower geometry is typical for wind pressures between 750 and 1500 N/m2 , with individual spans of 350 m, typical for 750 N/m2 wind pressure and 65 N/m2 of maximum allowed conductors tensile strength. Average ground flash density is taken at 1 (km−2 year−1 ) for the entire TL route. Tower geometry follows: h = 27 m, y = 24 m; span between tower consoles (arms) is 2 m; top console length a = 2.5 m, middle console length is 3 m and bottom console length is 3.5 m. Maximum phase conductor sag equals 3 m while that of the shield wire equals 2 m. Phase conductor DC resistance is 0.114 (/km) with 10.95 mm diameter, and that of the shield wire is 0.304 (/km) with 8 mm diameter. Insulator string (i.e. archorn) length equals 0.9 m. Fig. 2 has been obtained from (16) by means of applying the Armstrong and Whitehead EGM to the test case transmission line. On the left side of this figure is the quantile–quantile plot for the log-normal distribution and on the right side is the histogram of the variate generated from this distribution, with the PDF described by (16) superimposed on the histogram [43, Ch. 1]. The fitting of the distribution to the random data—generated from (16) by means of the Acceptance–Rejection method—is managed using the Maximum Likelihood Estimation of its parameters, e.g. [42, Ch. 2,43, Ch. 6]. This variate is further combined with a variate depicting the lightning current front duration, by means of the procedure described in Section 5, accounting at the same time for the correlation coefficient between these variates. Fig. 3 depicts a typical outcome of this procedure, in terms of the scatter plot with superimposed histograms for each of the variates, where the influence of the correlation between the variates is clearly evident [43, Ch. 3]. This has the effect of producing proportionately few lightning currents which at the same time have large amplitudes and short front durations. The lightning current parameters from the original set have been employed. Annually expected number of direct lightning strikes per 100 km of this TL is based on the expression (11), where the influence of the statistical distribution of lightning strikes along the TL span

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

253

Fig. 2. Log-normal distribution of lightning current amplitudes incident to the test case transmission line obtained using the Armstrong and Whitehead EGM (left: qq plot, right: histogram and PDF fit).

length has been taken into account by the coefficient Pt = 0.6. In Table 1 are presented estimated numbers of direct lightning strikes per 100 km-years of this TL, for several different electrogeometric models. Fig. 4 graphically depicts the BFOR probability (top plot) and BFOR per 100 km-years (bottom plot) as a function of tower

footing impedance, obtained for the test case TL with three different insulator flashover models, EGM according to Brown and Whitehead, and lightning current parameters from the original set. The tower footing impedance, used for the abscissa, is a low-frequency and low-current one, and can be obtained from the TL design documentation. The error bars on the top plot and “whiskers” on the

Fig. 3. Distribution of random variates from the original set of current amplitudes and front durations while accounting for their TL incidence and statistical correlation.

254

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

Table 1 Expected number of direct lightning strikes per 100 km-years of transmission line. Model

Strikes/100 km-years

Wagner Armstrong and Whitehead Brown and Whitehead Love Mousa and IEEE-1995

9.8 16.6 15.1 13.2 11.5

Fig. 5. The expected BFOR probability (top plot) and BFOR per 100 km-years (bottom plot) of the test case TL with three different EGM models and leader-progression model of the insulator flashover characteristic.

Fig. 4. The expected BFOR probability (top plot) and BFOR per 100 km-years (bottom plot) of the test case TL with three different insulator flashover models and lightningcurrent parameters drawn from the original set.

bars of the bottom plot of Fig. 4 come from the two-sided confidence interval of the BFOR probability, defined by (19), and determined for the 95% confidence level. It can be nicely observed from this figure that the leaderprogression model and equal-area model provide practically the same expected BFOR values, while V − t function model predicts higher values, which corroborates findings of other authors on similar transmission lines. In order to account for the influence of different EGMs on the transmission line BFOR, following three EGM models are used: Brown and Whitehead, Love and Mousa [5, Ch. 6]. Fig. 5 depicts the expected BFOR probability (top plot) and BFOR per 100 km-years (bottom plot) of the test case TL with three different EGM models applied and leader-progression model of the insulator flashover characteristic. In order to estimate the influence of the correlation coefficient (between the lightning current amplitudes and front durations) on the BFOR probability and BFOR per 100 km-years, following three values of the correlation coefficient—for the original set of lightning-current parameters—are selected: 0.1 (weak correlation), 0.47 (default correlation) and 0.9 (strong correlation). Additionally, the leader progression model is used for representing the insulator flashover and EGM according to Brown and Whitehead. Fig. 6 presents the expected BFOR probability (top plot) and BFOR per 100 km-years (bottom plot) for three different correlation levels between lightning current amplitudes and front durations. It is clearly evident from this figure that by neglecting the correlation one increases the BFOR probability. At the same time, strong correlation drastically decreases the BFOR probability. In order to assess the influence of the variation of statistical parameters of lightning currents incident to TL on the BFOR probability and BFOR per 100 km-years, the alternative set of lightning-current parameters is employed. Coefficient of

Fig. 6. The expected BFOR probability (top plot) and BFOR per 100 km-years (bottom plot) of the test case TL with three different correlation coefficients and lightningcurrent parameters drawn from the original set.

correlation was held at 0.47 level. The EGM according to Brown and Whitehead is applied. A typical distribution of random variates, depicting lightning current amplitudes and front durations from this alternative population, is graphically displayed in Fig. 7, in terms of the scatter plot with superimposed histograms for each of the variates [43, Ch. 3]. A drastic difference in TL incident lightning currents can be nicely observed by comparing Fig. 3 with Fig. 7. Now, there is proportionately far larger number of lightning currents having high amplitudes and short front durations. The BFOR probability and BFOR per 100 km-years, obtained from the alternative set of lightning-current parameters, is graphically depicted in Fig. 8. The influence of the lightning current parameters on the BFOR probability—all other parameters staying the same—is clearly evident and significant. This influence of the alternative set of lightning current parameters manifest itself through the increased expected BFOR on the test case TL, regardless of the insulator flashover model being applied, which can be nicely observed from Fig. 8.

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

255

Fig. 7. Alternative distribution of random variates depicting lightning current amplitudes and front durations while accounting for their TL incidence and statistical correlation.

Fig. 8. Test case transmission line BFOR probability (top plot) and BFOR per 100 kmyears (bottom plot) with alternative lightning data set and three different insulator flashover models.

The importance of the lightning current front duration on the BFOR probability can be further illustrated by holding it fixed on a certain value while randomly drawing amplitudes from the probability distribution depicted with (11), with I = 31.1 (kA) and  lnI = 0.484. For that purpose, Fig. 9 graphically displays the BFOR probability obtained from two different fixed values of lightning current front duration (with equal-area model of the insulator flashover and EGM according to Brown and Whitehead). The influence of the lightning-current front duration on the BFOR probability is quite evident from this depiction.

Fig. 9. Test case transmission line BFOR probability with two different fixed values of lightning current front duration.

The sensitivity study presented within this section reveals several factors influencing the expected BFOR on TLs and their relative importance. At the top of this list are the lightning-current parameters incident to transmission line. Here, a possibility of the statistical realisation of the combination of high lightning-current amplitudes with short front durations features prominently in increasing the expected BFOR, making the statistical correlation between them important for consideration. Next in line of importance is the model of the insulator string flashover, followed by the type of EGM of the TL lightning incidence being used. Tower surge impedance is of lesser importance in HV lines (due to their

256

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257

restricted height), although the BFOR generally increases with the increase of the tower height. It should be mentioned that the type of lightning current source used in simulations have influence on the BFOR; however, Heidler source is generally recommended in EMTP-ATP lightning surge transient analysis. It is noteworthy to mention that special lightning environment should be accounted for (e.g. notorious Japanese winter lightning) whenever possible in BFOR analysis, along with other TL features and keraunic levels. Particular importance is to be paid to the tower grounding system modelling, especially if it cannot be regarded as being concentrated (which has been invariably assumed throughout this paper), e.g. [16,22]. It is interesting to note that the discrepancy between simpler and more advanced models of the BFOR estimation can be as high as 100%, with a strong dependence on the tower grounding model being applied for the simulation [22]. Also, it could be noted here that several additional parameters influence the BFOR on TLs, not accounted-for here, such as: positive lightning strikes, subsequent strikes within the lightning flash, inclination of the lightning incidence path, upward answering leaders initiated by the descending step leader in lightning attachment process, effects of corona, etc. Furthermore, some aspects of the TL modelling for backflashover lightning transients are of concern and still under scrutiny: tower model, tower grounding model, insulator flashover model, and model of corona on TL conductors.

7. Conclusion This paper presented a comprehensive methodology for the analysis of the backflashover occurrence rates on HV transmission lines. Many of the major influential factors affecting the BFOR on HV transmission lines have been accounted for. It can be seen from the presented analysis that quite a few different parameters influence the backflashover occurrence on the HV transmission line towers, having various degrees of importance. The particular significance of the presented methodology could be seen in terms of the sensitivity analysis that it allows, which, in its application to the test case TL, provided levels of influences that different parameters exhibit in forming the backflashover occurrence rate. The downside of the presented methodology is in the long runtime needed for executing the successive EMTP-ATP simulations, necessitated by the application of the Monte Carlo method. The presented methodology, in general, could find its application in aiding the TL designer in optimising the tower geometry for different voltage levels and tower types. It can be of particular benefit in determining the probability of backflashover occurrence on the first several TL towers emanating from the HV substation, which is important in designing the substation overvoltage protection and selecting the parameters of HV metal-oxide surge arresters. Also, the presented methodology could find its application in determining the probability of backflashover occurrence on the several TL towers which are found to be exposed to frequent direct lightning strikes, e.g., part of the TL crossing a mountain ridge or traversing through a terrain associated with a high keraunic level. This analysis is important in selecting the TL towers where the surge arresters are to be installed for reducing the number of line outages due to backflashovers, i.e., optimising the number of surge arresters per tower, along with the total number of consecutive towers carrying the surge arresters. It could also be of the benefit in selecting the transmission line surge arrester parameters, particularly their energy absorption capability. Finally, employment of the methodology of this type (statistical simulation runs) could aid in improving the existing standards and recommendations in the field of TL shielding.

References [1] I.S. Grant, J.G. Anderson, A.R. Hileman, W. Janischewskyj, G.E. Lee, V.J. Longo, W. Chisholm, D. Parrish, N.K. Roukos, C.T. Main, E. Whitehead, J.T. Whitehead, A simplified method for estimating lightning performance of transmission lines, IEEE Trans. Power Apparatus Syst. 104 (4) (1985) 919–932, A report prepared by the Working Group on Lightning Performance of Transmission Lines, IEEE Power Engineering Society, Transmission and Distribution Committee. [2] J.T. Whitehead, W.A. Chisholm, IEEE working group report: estimating lightning performance of transmission lines II – updates to analytical models, IEEE Trans. Power Deliv. 8 (3) (1993) 1254–1267, A report of the Working Group on Estimating the Lightning Performance of Transmission Lines. [3] A. Eriksson, L. Dellera, J. Elovaara, A. Hileman, R. Houlgate, J. Huse, K.W. Weck, Guide to procedures for estimating the lightning performance of transmission lines, CIGRE Brochure 63 (1991) 1–61, Working Group 01 (Lightning) of Study Committee 33 (Overvoltages and Insulation Co-ordination). [4] C.A. Nucci, A survey on CIGRE and IEEE procedures for the estimation of the lightning performance of overhead transmission and distribution lines, in: Xth International Symposium on Lightning Protection, Curitiba, Brazil, 2009, pp. 151–165. [5] A.R. Hileman, Insulation Coordination for Power Systems, CRC Press, Boca Raton, FL, 1999. [6] H.W. Dommel, EMTP Theory Book, 2nd ed., Microtran Power System Analysis Corporation, Vancouver, BC, 1992. [7] EPRI, Electromagnetic Transients Program (EMTP): Application Guide (Tech. rep.), Electric Power Research Institute, EMTP Development Coordination Group, EPRI EL-4650, 1988. [8] A. Ametani, T. Kawamura, A method of a lightning surge analysis recommended in Japan using EMTP, IEEE Trans. Power Deliv. 20 (2) (2005) 867–875. [9] IEEE WG, Modeling guidelines for fast front transients, IEEE Trans. Power Deliv. 11 (1) (1996) 493–506, Report Prepared by the Fast Front Transients Task Force of the IEEE Modeling and Analysis of System Transients Working Group. [10] CIGRE WG, Guidelines for representation of network elements when calculating transients, CIGRE Brochure 39 (1988) 1–29, Working Group 33.02 (Internal Overvoltages) of Study Committee 33 (Overvoltages and Insulation Co-ordination). [11] J.A. Martinez-Velasco (Ed.), Power System Transients: Parameter Determination, CRC Press, Boca Raton, FL, 2010. [12] J.A. Martinez, F. Castro-Aranda, Lightning performance analysis of overhead transmission lines using the EMTP, IEEE Trans. Power Deliv. 20 (3) (2005) 2200–2210. [13] M.Z.A. Ab Kadir, I. Cotton, Application of the insulator coordination gap models and the effect of line design to backflashover studies, Electr. Power Energy Syst. 32 (2010) 443–449. [14] T. Ito, T. Ueda, H. Watanabe, T. Funabashi, A. Ametani, Lightning flashovers on 77-kV systems: observed voltage bias effects and analysis, IEEE Trans. Power Deliv. 18 (2) (2003) 545–550. [15] F.M. Gatta, A. Geri, S. Lauria, Backflashover simulation of HV transmission lines with concentrated tower grounding, Electr. Power Syst. Res. 73 (2005) 373–381. [16] F.M. Gatta, A. Geri, S. Lauria, M. Maccioni, Backflashover simulation of HV transmission lines with enhanced counterpoise groundings, Electr. Power Syst. Res. 79 (2009) 1076–1084. [17] W. Nowak, R. Tarko, Computer modelling and analysis of lightning surges in HV substations due to shielding failure, IEEE Trans. Power Deliv. 25 (2) (2010) 1138–1145. [18] I.F. Gonos, L. Ekonomou, F.V. Topalis, I.A. Stathopulos, Probability of backflashover in transmission lines due to lightning strokes using Monte-Carlo simulation, Electr. Power Energy Syst. 25 (2003) 107–111. [19] R. Holt, T.T. Nguyen, Monte Carlo estimation of the rates of lightning strikes on power lines, Electr. Power Syst. Res. 49 (1999) 201–210. [20] S.J. Shelemy, D.R. Swatek, Monte Carlo simulation of lightning strikes to the Nelson River HVDC transmission lines, in: International Conference on Power System Transients, Rio de Janeiro, Brazil, 2001. [21] A. Shafaei, A. Gholami, R. Shariatinasab, Probabilistic evaluation of lightning performance of overhead transmission lines, considering non-vertical strokes, Sci. Iran. D 19 (3) (2012) 812–819. [22] F.M. Gatta, A. Geri, S. Lauria, M. Maccioni, A. Santarpia, An ATP-EMTP Monte Carlo procedure for backflashover rate evaluation: a comparison with the CIGRE method, Electr. Power Syst. Res. 113 (2014) 134–140. [23] A. Ametani, Lightning surge analysis by EMTP and numerical electromagnetic analysis method, in: 30th International Conference on Lightning Protection, Cagliary, Italy, 2010, pp. 1–18, Invited Lecture L1. [24] J.A. Martinez, B. Gustavsen, D. Durbak, Parameter determination for modeling system transients—Part I: Overhead lines, IEEE Trans. Power Deliv. 20 (3) (2005) 2038–2044, IEEE PES Task Force on Data For Modeling System Transients of IEEE PES Working Group on Modeling and Analysis of System Transients Using Digital Simulation. [25] L. Prikler, H.K. Hoidalen, ATPDRAW version 5.6 for Windows 9x/NT/2000/XP/Vista: User’s Manual, Norwegian University of Technology, Trondheim, Norway, 2009. [26] M.A. Sargent, M. Darweniza, Tower surge impedance, IEEE Trans. Power Apparatus Syst. 88 (5) (1969) 680–687. [27] T. Suzuki, K. Miyake, Experimental study of breakdown voltage-time characteristics of large air gaps with lightning impulses, IEEE Trans. Power Apparatus Syst. 96 (1) (1977) 227–233.

P. Sarajcev / Electric Power Systems Research 119 (2015) 247–257 [28] T. Shindo, T. Suzuki, A new calculation method of breakdown voltage-time characteristics of large air gaps, IEEE Trans. Power Apparatus Syst. 104 (6) (1985) 1556–1563. [29] H. Motoyama, Experimental study and analysis of breakdown characteristics of long air gaps with short tail impulse, IEEE Trans. Power Deliv. 11 (2) (1996) 972–979. [30] M. Kizilcay, C. Neumann, Backflashover analysis for 110-kV lines at MultiCircuit overhead line towers, in: International Conference on Power System Transients, Lyon, France, Paper No. 293, 2007. [31] R. Anderson, A. Eriksson, Lightning parameters for engineering application, Electra 69 (1980) 65–102. [32] CIGRE, Lightning Parameters for Engineering Applications (Tech. rep.), CIGRE, Working Group C4.407, 2013. [33] IEEE WG, Parameters of lightning strokes: a review, IEEE Trans. Power Deliv. 20 (1) (2005) 346–358, IEEE Working Group of the Lightning and Insulator Subcommittee of the Transmission and Distribution Committee. [34] A. Borghetti, C.A. Nucci, M. Paolone, Estimation of the statistical distributions of lightning current parameters at ground level from the data recorded by instrumented towers, IEEE Trans. Power Deliv. 19 (2004) 1400–1409. [35] J. Takami, S. Okabe, Observational results of lightning current on transmission towers, IEEE Trans. Power Deliv. 22 (1) (2007) 547–556. [36] V. Cooray (Ed.), Lightning Protection, IET Power and Energy Series, vol. 58, The Institution of Engineering and Technology, London, UK, 2010.

257

[37] V. Cooray, Lightning interception: a review of simulation procedures utilized to study the attachment of lightning flashes to grounded structures, Electra 257 (2011) 48–55, on behalf of CIGRE Working Group C4.405. [38] P.N. Mikropoulos, T.E. Tsovilis, Lightning attachment models and maximum shielding failure current of overhead transmission lines: implications in insulation coordination of substations, IET Gener. Transm. Distrib. 4 (12) (2010) 1299–1313. [39] R. Lambert, A. Xemard, G. Fleury, E. Tarasiewicz, A. Morched, Probability density function of the lightning crest current at ground level: estimation of the lightning strike incidence on transmission line, in: International Conference on Power System Transients, Budapest, Hungary, 1999. [40] L. Ekonomou, I.F. Gonos, I.A. Stathopulos, Lightning performance assessment of Hellenic high voltage transmission lines, Electr. Power Syst. Res. 78 (2008) 703–712. [41] A.M. Mousa, K.D. Srivastava, The distribution of lightning strokes to towers and along the span of shielded and unshielded power lines, Can. J. Electr. Comput. Eng. 15 (3) (1990) 115–122. [42] W. Hauschild, W. Mosch, Statistical Techniques for High-Voltage Engineering IEE Power Series, vol. 13, Peter Peregrinus Ltd, London, UK, 1992. [43] W.K. Hardle, L. Simar, Applied Multivariate Statistical Analysis, 3rd ed., Springer-Verlag, Berlin, Germany, 2012.