EMTP IMPLEMENTATION OF A MONTE CARLO METHOD FOR LIGHTNING PERFORMANCE ANALYSIS OF TRANSMISSION LINES

An accurate calculation of lightning overvoltages is an important issue for the analysis and design of overhead transmission lines. The different parts of a transmission line that are involved in lightning calculations must be represented taking into account the frequency ranges of transients associated to lightning. In addition, the procedures to be used in these calculations must be developed considering the random nature of lightning phenomena. Several simulation tools have been used to estimate the lightning performance of transmission lines. The most popular approaches are those based on a time-domain simulation technique for which adequate procedures and transmission line models have to be developed. This paper presents a summary of the computational efforts made by the authors for the development and implementation in an EMTP-like tool of a Monte Carlo procedure, as well as the models of some transmission line components, aimed at analyzing the lightning performance of transmission lines. An actual test line is used to illustrate the scope of this procedure and the type of studies that can be performed.


INTRODUCTION
Lightning is a physical phenomenon of random nature, and one of the main causes of overhead line failures.The lightning performance of an overhead line can be measured by a ashover rate, usually expressed as the number of ashovers per 100 km per year [1,2].Transmission lines are usually shielded by one or several wires; therefore, lightning failures can be caused by strokes to either a shield wire or a phase conductor.Overvoltages induced by strokes to ground can be neglected.Shielding failures cannot be totally prevented, but the number of strokes to phase conductors is usually very low.The ashover rate of a transmission line is therefore divided into the Back ashover Rate (BFOR) and the Shielding Failure Flashover Rate (SFFOR).To obtain both quantities an incidence model is required to discriminate strokes to shield wires from those to phase conductors and those to ground.
Due to the random nature of lightning, the ashover rate calculation must be based on a statistical approach.On the other hand, lightning overvoltage calculations are generally performed by applying a time-domain technique and using adequate models for the frequency ranges associated to lightning transients, although other approaches based on a frequency-domain technique can be also considered [3][4][5].
The application of a Monte Carlo method is the usual procedure method for this purpose.Such procedure may consist of the following steps: generation of random numbers to obtain those parameters of the lightning stroke and the overhead line of random nature; application of an incidence model to deduce the point of impact of every lightning stroke; calculation of the overvoltage generated by each stroke, depending on the point of impact; calculation of the ashover rate [6].This paper is aimed at detailing the main computational aspects of the procedure and models developed by the authors for assessment of the lightning performance of overhead transmission lines [6][7][8][9].A Monte Carlo procedure and new models for representing several parts of a transmission line in lightning overvoltage calculations have been implemented in the ATP package.ATP is an acronym that stands for Alternative Transients Program, and it is a well known member of the EMTP (ElectroMagnetic Transients Program) family [10].
This document has been organized as follows.The next section includes a summary on the models used for representing the transmission line in lightning overvoltage calculations.The four subsequent sections present respectively lightning characterization, as assu-med in this paper, the main features of the developed Monte Carlo procedure, ATP solution methods and the capabilities of the procedure implemented in ATP.The results derived from the application of the Monte Carlo procedure as well as the main conclusions of this paper are presented in the last sections.

MODELING FOR LIGHTNING OVERVOLTAGE CALCULATIONS
The models used to represent the different parts of the transmission line involved in overvoltage calculations are detailed in the following paragraphs.For more information on the guidelines to be considered for representation of overhead transmission lines in lightning studies see [11][12][13].
1) Shield wires and phase conductors of the transmission line can be modeled by three or more spans at each side of the point of impact.A rigorous representation of each span should be based on a multi-phase frequencydependent untransposed distributed-parameter line model.However, for lightning over-voltage calculations, a constant parameter line model can be accurate enough [12].In this paper line parameters are constant and calculated at 500 kHz.2) The line termination at each side of the above model, needed to avoid re ections that could affect the simulated overvoltages around the point of impact, is represented by a long enough section, whose parameters are also calculated at 500 kHz.3) Several models have been proposed to represent transmission line towers; they have been developed using a theoretical approach or based on an experimental work [7].The simplest representation is a lossless distributed-parameter transmission line, characterized by a surge impedance and a travel time.This model can suf ce for lines with towers shorter than 50 m [1].4) Phase voltages at the instant at which the lightning stroke impacts the line have to be included in calculations; their values are deduced by randomly determining a phase reference angle.5) A lightning stroke is represented as an ideal current source (in nite parallel impedance) whose parameters, as well as its polarity and multiplicity, are randomly determined according to the distribution density functions recommended in the literature [1,[14][15].See next section for more details.6) The representation of insulator strings relies on the application of the leader progression model [1,[16][17], which can be used to account for non-standard lightning voltages.According to this model, the ashover mechanism consists of three steps: corona inception, streamer propagation and leader propagation.When the applied voltage exceeds the corona inception voltage, streamers propagate along the insulator string; if the voltage remains high enough, these streamers will become a leader channel.A ashover occurs when the leader crosses the gap between the cross-arm and the conductor.Therefore, the total time to ashover can be expressed as follows where t c is the corona inception time, t s is the streamer propagation time and t l is the leader propagation time.Usually t c is neglected because it is very short in comparison to the other two times.According to [17], t s can be calculated as follows where E 50 is the average gradient at the critical ashover voltage and E is the maximum gradient in the gap before breakdown.The leader propagation time, t l , can be obtained from the following equation where V(t) is the voltage across the gap, g is the gap length, l is the leader length, E l0 is the critical leader inception gradient, and k l is a leader coef cient.The leader propagation stops if the gradient in the unbridged part of the gap falls below E l0 .
7) Footing impedance modeling is a critical aspect.
A nonlinear frequency-dependent representation is required to obtain an accurate simulation [18,19].Since the information needed to derive such a model is not always available, a lumped nonlinear resistance is usually chosen for representing the footing impedance, although it cannot be always adequate.The value of this resistance is approximated by the following expression [1,11,20] where R o is the footing resistance at low current and low frequency, I is the stroke current through the resistance, and I g is the limiting current to initiate suf cient soil ionization, which is given by being the soil resistivity (ohm-m) and E 0 the soil ionization gradient (400 kV/m, [19]).

LIGHTNING STROKE CHARACTERIZATION
The following paragraphs detail the most important aspects to be considered for a full characterization of lightning ashes.For a more complete characterization and an updated bibliography on this subject see [14] and [15].

Polarity
Most lightning ashes are of negative polarity.The incidence of positive ashes increases during winter, although very rarely their percentage exceeds 10% [14,15].

Multiplicity
Negative ashes can consists of multiple strokes, while positive ashes have usually a single stroke.Very different percentages of single-stroke ashes have been reported in the literature, see [14,[21][22].Less than 10% of positive ashes have multiple pulses [15]; however, since the number of recorded positive ashes is too low, it is usually assumed they are single-stroke ashes.

Return-Stroke Waveform
A concave waveform, with no discontinuity at t=0, is an accurate representation of the wave front of a negative return stroke.Several expressions have been proposed for such a form, being the so-called Heidler model one of the most widely used, see gure 1 [23].It is given by the following expression: where I p is the peak current, is a correction factor of the peak current, n is the current steepness factor, k = t/ 1 , while 1 , 2 are respectively time constants determining current rise and decay time [23].Ingeniare.Revista chilena de ingeniería, vol.16 Nº 1, 2008 Main parameters used to de ne this waveform in the present work are the peak current magnitude, I 100 , the rise time, t f (= 1.67 (t 90 -t 30 )), and the tail time, t h , i.e. the time interval between the start of the wave and the 50% of the peak current on tail.
Reported waveforms for the rst negative stroke show some differences in the peak zone with respect the waveform depicted in gure 1 [15,16].Negative subsequent strokes do not show a very pronounced concavity of the wavefront [15], and they can be represented with enough accuracy by the Heidler model.Although positive return strokes do not have enough common features to produce an acceptable waveform [15], the above waveform can be also used to represent them.

Probability Distribution of Return-Stroke Parameters
The statistical variation of the lightning stroke parameters can be approximated by a log-normal distribution, with the following probability density function [15] where lnx is the standard deviation of lnx, and x m is the median value of x.
The joint probability density function of two stroke parameters can be expressed as follows exp .If x and y are independently distributed, then c =0, and p(x,y)=p(x)p(y).

Return-Stroke Parameters
Although negative ashes have multiple strokes, only the rst and the second strokes are of concern for transmission insulation levels; i.e. peak current magnitudes of the third and the subsequent strokes are much lower than those of the previous strokes and they do not represent a threat to transmission lines.In addition, only an accurate knowledge of the parameters of the rst negative stroke can be crucial for transmission lines with rated voltage 400 kV and above [8].Since the test line analyzed in this paper is an actual 400 kV transmission lines, only results derived from assuming single-stroke ashes are presented.
The peak current magnitudes of positive strokes are larger than those of negative polarity, while their front and tail times are much longer.In addition, they present a seasonal variation: the number of positive ashes increases during winter [24,25], when their statistical parameters are very different from those of negative ashes.However, the same conclusion is derived when comparing winter positive flashes to summer positive flashes, whose parameters are similar to those of negative strokes.Due to this seasonal variation and to the uncertainties related to the percentages of positive and negative ashes, the calculation of flashover rates will be performed by separating positive and negative polarity ashes; that is, only positive or negative polarity ashes will be assumed when estimating ashover rates.

MONTE CARLO PROCEDURE
The following paragraphs detail the most important aspects of the procedure implemented in ATP, for the analysis of the lightning performance of transmission lines [6].
1) The values of random variables and parameters are estimated.This step includes the calculation of the parameters of the lightning stroke (peak current, rise time, tail time, and location of the leader channel), the phase conductor voltages, the footing resistance and the insulator strength.
2) The last step of a return stroke is determined by the electrogeometric model, see gure 2 [1].The values of the striking distances are calculated according to the following expressions where and are constants that depend on the object, is a constant that depend on the electrogeometric model and I p is the stroke current [1,2].3) Overvoltage calculations are performed once the point of impact has been determined.The only difference between models for back ash and shielding failure simulations is the node to which the current source that represents the stroke must be connected.Overvoltages caused by nearby strokes to ground are not simulated, since their effect can be neglected for transmission insulation levels.4) If a ashover occurs in an insulator string, the counter is increased and the ashover rate updated.5) The convergence of the Monte Carlo method is checked; the procedure is stopped when the calculated probability density functions of all random variables match their theoretical functions within the speci ed error.

ATP SOLUTION METHODS, CAPABILITIES AND APPLICATIONS
The ATP was originally developed for simulation of electromagnetic transients in power systems.A transient simulation is carried out with a xed time step selected by the user, using the trapezoidal rule of integration and the Bergeron's method.However, the program can also be used to perform AC steady state calculations, to obtain system impedance as a function of the frequency and to calculate harmonic power ows.Solution methods to solve systems with nonlinear components have been also implemented [10].
The program can represent control systems and interface them with an electric network.Several non-simulation supporting routines are also available to create model data les.These supporting routines can be used for computation of line parameters or derivation of coupled RL matrix aimed at representing multi-phase, multiwinding transformers.
The ATP package is integrated by at least three tools: ATPDraw, a graphical preprocessor for creating/editing input les; TPBIG, the main engine for transients and harmonics simulations; one postprocessor for plotting simulation results.ATPDraw is an interactive program that can act as a shell for the whole package: users can control the execution of all programs integrated in the package from ATPDraw.Standard components are supported by this tool, although users can create custommade objects based on ATP capabilities.Figure 3 shows a schematic diagram of the connectivity between simulation capabilities, supporting routines, external programs and all types of les.
ATP capabilities allows users to simulate a data case several times before deactivating the program, change parameters of the system under simulation according to a given law, disconnect or activate some components, and carry out some calculations by using external programs.In addition, it is possible to modify "on line", if required, the simulation time or the number of runs.
The usage of the ATP package is based on the following concepts [26]: 1) Multiple run option: A data case can be simulated as many times as needed, while changes are introduced into the system at every run.This option is known as Pocket Calculator Varies Parameters (PCVP) [27].
It can be used to perform statistical and sensitivity studies.However, it can also be used in many other applications.After the target of the study has been set, this option can be used to run the case as much as required while one or several parameters are gradually adjusted or the system topology is modi ed, until the target is reached.2) Open system: A link to external programs/tools can be established before, during and after a simulation to take advantage of the capabilities of these tools and to add or test new capabilities.This option can be used, for instance, to link the ATP to MATLAB and take advantage of its toolboxes, or to run a custom-made program that can derive the parameters of a power component using a data conversion procedure not yet implemented in the package.3) Data symbol replacement: $PARAMETER is a declaration that can be used to replace data symbols of arbitrary length prior to a simulation [27].Up to three replacement modes can be used: simple character replacement (one string is replaced by other with the same length), mathematical replacement (string is replaced by a number deduced from a mathematical formula), integer serialization (used to encode strings within a DO loop).Conditional branching (IF-THEN-ENDIF) is a built-in feature that can be used to select between two or more choices.4) Data module: ATP-coded templates have been used in the past for the development of data modules that could facilitate the tool use by beginners, or to simplify the use of power components and extend modelling capabilities to more complex equipment [27].Custommade models are represented by a module and its associated ATPDraw icon.The development of new modules relies on the routine Data Base Module, but other ATP capabilities can be used to perform simple calculations with module arguments, to decide what parts of a module can be activated at a given run, or what parts should be deactivated.5) Custom-made model/routine: They are based on a feature known as TACS Type-69.Although it was developed for implementation of compiled custommade functions that could be embedded into a control unit, it can also be used for the development of new routines (e.g. the calculation of random variables) or models (e.g. an insulator string representation based on the leader propagation model).6) Interactivity: Several simulation modules will be usually involved in a general procedure.Interactivity between them is critical as calculations will be performed in several modules.The connectivity between a power system and a control section to pass variables in both ways has been an ATP feature since the earliest development of control capabilities.However, it has been the possibility of passing also parameters what has added exibility to some of the capabilities described above and increased the type of applications.The applications that can be covered by the ATP can be grouped as follows [26]: Time-domain simulations: They are generally used for simulation of transients, such as switching or lightning overvoltages.However, they can also be used for analyzing harmonic distortion.
Frequency-domain simulations: ATP capabilities can be used to obtain the driving point impedance at a particular node versus frequency, detect resonance conditions, design lter banks or analyze harmonic propagation.
Sensitivity analysis: It is usually performed to evaluate the relationship between variables and parameters.
When one or more parameters cannot be accurately speci ed, this analysis will determine the range of values for which they are of concern.
Statistical studies: Several ATP capabilities can be combined to perform all types of Monte Carlo simulations.

IMPLEMENTATION OF THE MONTE CARLO PROCEDURE
Figure 4 shows the owchart of the procedure implemented in ATP when the lightning performance of the test line is analyzed by assuming that only single-polarity singlestrokes are produced.The ATP capabilities that were used in the development and application of this procedure are listed below.
The multiple-run option is used to perform all the runs required by the Monte Carlo method.The procedure can be stopped when either the convergence is achieved or the maximum number of runs is reached.A compiled routine was developed and embedded into the source code to obtain the values of the random variables/parameters that must be generated at every run according to the probability distribution function assumed for each one.The generation of random numbers is based on a linear congruential multiplicative generator with a loop of sequences greater than 2 144 [28,29].Table 1 lists the type of function assumed for each variable or parameter of random nature.For more information on the algorithms to be used for each density function, see [30].Doing so, the routine can be treated as a built-in capability.
Phase-to-ground voltages across insulators are continuously monitored; when the voltage stress in a single phase exceeds the strength, the ashover counter is increased and the simulation is stopped.That is, every run has a different simulation time, depending on insulator string behavior and the maximum front time of lightning strokes.

Footing resistance Normal
The Monte Carlo procedure is stopped when the maximum number of runs, speci ed by the user, is reached or when the probability density functions of all the random variables match their theoretical functions within the speci ed error.In this work, the resulting and theoretical distributions are compared at 10, 30, 50, 70 and 90% of the cumulative distribution functions.More than 10000 runs were needed to match them within an error margin of 10%.For an error margin of 5% no less than 30000 runs were needed (see next section).
A report showing the main input and output variables, as well as the progress of the ashover rate, can be printed with a frequency chosen by the user.
Figure 5 shows the list of compiled routines and modules developed for application of these routines.The interactions between the different routines and tasks are depicted in gure 6.
Note that there are two types of routines: those developed for application of the Monte Carlo method (generation of random numbers, estimation of different random variables and parameters, application of the electrogeometric model, control of the convergence and information output) and those developed for representation of some parts of a transmission line (insulator and footing resistance models), that were not available as built-in models.
A module or library element has been created for representing a transmission line tower using the ATP capability known as DATA BASE MODULE.The module is stored as an ordinary le and used every time the component has to be simulated using an $INCLUDE command, with a list of arguments describing the tower geometry.A component (e.g. the tower) is represented by a template that consists of different network elements (e.g.branches, switches or sources) that are needed to build the equivalent circuit of that component.The same capability has been used to create the modules listed in the left column of gure 5; in this case the modules work as an interface that facilitates the usage of the new built-in routines and models.Figure 6.Interactions between tasks implemented for lightning ashover rate calculation.

Test Line
Figure 7 shows the tower design for the line tested in this paper.It is a 400 kV line, with two conductors per phase and two shield wires, whose characteristics are provided in table 2.

Transmission Line and Lightning Parameters
A model of the test line was created using ATP capabilities and following the guidelines summarized in the section on modeling guidelines.
The line was represented by three 400-m spans plus a 30-km section as line termination at each side of the point of impact.
The surge impedance of towers was calculated according to the expression suggested by CIGRE [17].A value of 134 was estimated for all towers.
Table 3 shows the parameters to be used in the insulator equations that were detailed above.The value of the average gradient at the critical ashover voltage, E 50 , was assumed to be the same that E l0 for both negative and positive ashes.The striking distance of insulator strings was in all cases 3.212 m.
Only single stroke ashes were considered.A return stroke was represented by a concave waveform, with parameter n, to be speci ed in equation ( 6), equal to 5.
The following probability density functions were assumed for each random value: Lightning stroke: Peak current magnitude, rise time and tail time were determined by assuming a log-normal distribution for all of them.Table IV shows the values used for each parameter.In addition, parameters were independently distributed, i.e. c =0.Since unrealistic values of the peak current magnitude can be generated using these parameters, the distributions have been truncated for both polarities at 500 kA.No ashovers other than those across insulator strings have been considered.

Simulation Results
The in uence of the ash polarity was analyzed by assuming that all ashes were either only negative or only positive.After 40000 runs with each polarity, the ashover rate was respectively 0,6025 and 0,9075 per 100 km-year for negative and positive ashes [8].These results were deduced by assuming that the ash ground density was N g = 1 ash per km 2 and year.
Figure 8 shows the results obtained with negative polarity ashes.Similar results, with different range of values, were obtained with positive polarity ashes.From these results the following conclusions can be derived: There is a range of peak current magnitudes for each type of ashover and a range of values that do not cause ashover.Since the parameters assumed for each polarity were different, the distributions are also different, although back ashovers are not caused with peak current magnitudes below 100 kA, neither with negative nor with positive strokes.

Discussion
The value of the rise time with which back-ashovers are caused shows a different range for each polarity: no back ashovers were caused by negative strokes with a rise time above 5 s, while there were back ashovers caused by positive strokes with a rise time longer than 10 s.The main factors that contribute to the latter result are the insulator strength, which is lower under positive strokes, and the tail duration of these strokes, which is much longer than that of negative strokes.According to the leader progression model, a ashover can be caused after the peak voltage if the tail is long enough.
Several parts of the transmission line model are not accurate enough: corona effect is not included in the line span models, voltages induced by the electric and magnetic elds of lightning channels to shield wires and phase conductors are neglected, the footing impedance model is too simple.In addition, a vertical direction is assumed for the stroke leader when it approaches earth.However, the incorporation of these aspects will not signi cantly vary the ashover rate, see for instance [9].In addition, it is worth mentioning that no study with a more realistic representation of all these aspects has been performed to date.Ingeniare.Revista chilena de ingeniería, vol.16 Nº 1, 2008 Sensitivity studies can be very useful to analyze the in uence of line and stroke parameters, and to determine what range of values might be of concern.
Although the number of parameters involved in lightning calculations is very high, only some of them can be accurately speci ed from the line geometry.For results derived from sensitivity studies presented in previous papers see [6][7][8][9].From those studies it is obvious that the ashover rate increases with the peak current magnitude and decreases with the rise time; in addition, the in uence of the tower footing resistance is not critical for low values of the soil resistivity and low values of R m .
A non-zero correlation coef cient between the probability density functions of the peak current magnitude and the rise time has been suggested [17].Figure 9 shows the relationship between ashover rate and correlation coef cient for negative ashes.A similar trend is deduced for positive ashes.According to results previously presented, the BFOR is very sensitive to the coef cient of correlation between the peak current magnitude and the rise time, while the SFFOR remains practically constant, irrespective of the value of c ; consequently, the Total Flashover Rate decreases as the coef cient of correlation increases.These results con rm that an accurate knowledge of this parameter is an important issue for ashover rate calculations.

CONCLUSIONS
EMTP-like programs are powerful and exible tools that can be successfully applied to the lightning performance analysis of overhead lines.This paper has shown how the ATP can be adapted for the study of transmission lines.
ATP capabilities have been applied to the development and implementation of a Monte Carlo procedure and several models that are adequate for the lightning performance analysis of transmission lines.The paper has presented the ATP features that can be used for such purpose, summarized the models needed for representing the different parts of a transmission line, shown how the ATP features can be applied to the development of some of them, illustrated the scope of the procedure and the models developed for the present work.
Although the transmission line model used in this work has some limitations, it is in the characterization of lightning where the limitations are more important.In fact, most EMTP-like tools have capabilities that allow users to develop models more accurate and sophisticated than some of those used in this work; however, the uncertainties with which lightning strokes are represented and the fact that modern detection systems cannot provide good estimates of all lightning parameters are more obvious and important limitations [31][32][33].

Figure 1 .
Figure 1.Parameters of a return stroke -Concave waveform.
c the coef cient of correlation.

Figure 2 .
Figure 2. Application of the electrogeometric model.

Figure 3 .
Figure 3. ATP simulation modules and supporting routines.

Figure 4 .
Figure 4. Diagram of the Monte Carlo procedure.

Figure 5 .
Figure 5. Models and routines implemented into the ATP.
current magnitude of strokes to shield wires that caused ashover.current magnitude of strokes to phase conductors that caused ashover.

Table 4 .
Statistical parameters of return strokes.