Lightning performance analysis of transmission lines using the Monte Carlo method and parallel computing

An accurate calculation of lightning overvoltages is crucial for the design of overhead lines. Given the random nature of the lightning phenomenon, the procedure to be used must be statistical and generally based on the Monte Carlo method. Several simulation tools have been used to date for estimating the lightning performance of transmission lines. Some popular approaches use an EMTP-like tool. This type of approach has some drawbacks: the line must be represented using a complex and sophisticated model, and the application of a Monte Carlo method requires a high number of runs. A solution to this problem is the usage of a parallel computing environment, which can reduce the simulation time in a ratio close to the number of cores used in calculations. This paper presents a summary of the effort made by the authors for the development of a Monte Carlo procedure aimed at estimating the lightning performance of transmission lines by means of an EMTP-like tool and the use of a multicore installation, in which the EMTP is controlled from a custom-made application implemented in MATLAB.


INTRODUCTION
Lightning is one of the main causes of overhead line failures.The lightning performance of an overhead line is usually measured by the number of flashovers per 100 km and year [1].One or two wires usually shield transmission lines.Since the level of overvoltages induced by strokes to ground is too low for transmission-level overhead lines, lightning failures can be due to strokes to either a shield wire or a phase conductor.Shielding failures cannot be totally prevented, but the number of strokes to phase conductors is usually very low.The lightning flashover rate (LFOR) of a transmission line is estimated by adding the Backflashover Rate (BFOR) and the Shielding Failure Flashover Rate (SFFOR) [1][2].To obtain both quantities an incidence model is required to discriminate strokes to shield wires from those to phase conductors and those to ground [1][2].
Due to the random nature of lightning, the flashover rate calculation must be based on a statistical approach.On the other hand, an accurate estimation of lightning overvoltages can be obtained by applying a time-domain technique and using adequate models for the frequency ranges associated to lightning transients.The application of a Monte Carlo method is the usual solution for this type of studies [3].However, the Monte Carlo method has a clear disadvantage: a high number of runs/ simulations is required to achieve accurate results.Besides, this is aggravated by the fact that the line model will be very complex and sophisticated [4].An obvious solution to this disadvantage is the usage of a multicore installation.The distribution of the Monte Carlo runs between several cores can reduce the computing time in a ratio close to the number of cores.In other words, the procedure can be fast if the number of cores is high.This paper details the main aspects of a procedure implemented for assessing the lightning performance of overhead transmission lines when the line model is implemented in the ATP (Alternative Transients Program) version of the EMTP [5].The Monte Carlo procedure has been implemented in a MATLAB application that pre-calculates random values, distributes the ATP tasks between the various cores involved in lightning overvoltage calculations, and post-process simulation results.This document has been organized as follows.The next two sections summarize the characterization of lightning strokes and the modeling guidelines used for representing transmission lines in lightning overvoltage calculations.Since these aspects have been already covered, see references [4] and [6], only a summary is provided here.The main contributions of this work are presented in two sections dedicated to respectively detail the implementation of the parallel MATLAB-ATP procedure and its application to the analysis of an actual overhead transmission line.The main conclusions are summarized in the last section.

LIGHTNING STROKE CHARACTERIZATION
The most important aspects for a full characterization of a lightning flash are summarized below.For more details see [7].
Polarity: Most lightning flashes are of negative polarity.The incidence of positive flashes increases during winter, although very rarely their percentage exceeds 10% [7].
Multiplicity: Negative flashes can consist of multiple strokes, while positive flashes have usually a single stroke, see references [7] and [8].Less than 10% of positive flashes have multiple pulses.
Return stroke waveform: A concave waveform, with no discontinuity at t = 0, is an accurate representation 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 [9].The waveform selected in this paper for representing the stroke current is a concave waveform based on a double exponential, see Figure 1, and defined by the following expression: ( ) with where I p is the peak current, t s is the start time of the waveform, t w = t when i 1 (t) = 0.999I p , and The parameters used to define this waveform (see Figure 1) are the peak current magnitude, I 100 = I p , the rise time, t f , and the tail time, t h (i.e. the time interval between the start of the wave and the 50% of the peak current measured on tail).Parameter n is a correction factor of the wave front.In this work, parameter n is fitted to fulfill t f = 1.67 (t 90 -t 30 ), being t 90 and t 30 the instants at which the stroke current reach respectively the 90% and the 30% of the peak current value; the resulting value is n = 3.24.The waveform shown in Figure 1 is based on that proposed in standards and used by some authors; see, for instance, [1].
Return stroke parameters: Although negative flashes have multiple strokes, only the first and the second strokes are of concern for transmission insulation levels; that is, 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.Actually, only an accurate knowledge of the parameters of the first negative stroke can be crucial for transmission lines with rated voltage 400 kV and above [10].
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 exhibit a seasonal variation: the number of positive flashes increases during winter [11][12], when their statistical parameters are very different from those of negative flashes.The same conclusion is derived when comparing winter positive flashes to summer positive flashes, whose parameters are similar to those of negative strokes.
The statistical variation of the lightning stroke parameters can be approximated by a log-normal distribution, with the following probability density function [7]: where s lnx is the standard deviation of lnx, and x m is the median value of x.
A non-zero correlation coefficient between the probability density functions of the peak current magnitude and the rise time is usually accepted [7].See also references [3] and [6].

MODELING FOR LIGHTNING OVERVOLTAGE CALCULATIONS
Modeling guidelines for representing the different parts of the transmission line involved in lightning overvoltage calculations are summarized in the following paragraphs.For more information on the guidelines see [4], [13][14][15].
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 frequency-dependent untransposed distributed-parameter line model.For lightning overvoltage calculations, a constant parameter line model with parameters calculated at a high frequency (e.g.500 kHz) can be accurate enough [4].2. The line termination at each side of the above model can be represented by a long-enough line section to avoid during the simulation time reflections that could affect the calculated overvoltages around the point of impact.Note that only waves reflected at the towers closer to the point of impact should affect the voltages caused by the lightning stroke current.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 [16].The simplest representation is a lossless distributedparameter transmission line, characterized by a surge impedance and a travel time.For an introduction to tower modeling see [17].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. Corona effect can impact the propagation of overvoltage surges associated with lightning strokes: corona introduces a time delay to the wave front [18].This time delay takes effect above the corona inception voltage and varies with surge magnitude.This variation can be expressed as a voltage-dependent capacitance added to the geometrical capaci-tance of the transmission line.In general, corona does not affect the tail of the surge.
Corona models, such as that shown in Figure 2, can be used to model the dynamic capacitance region of the q/V curve in a piecewise linear fashion.However, this model has some limitations: • It is based on lumped elements and must be lumped at sufficiently small intervals along the line to minimize the error introduced by the discretization.A minimum interval length of 50 meters has been suggested; however, shorter intervals can be advisable.• The model does not adequately address corona in a multiphase model.Here, the voltage dependence of the corona should be transformed into charge dependence, because corona depends on the electric field around the conductor.Corona effect is independent of the conductor size and geometry for high magnitude surges.For low magnitude surges, the effect will depend on conductor sizes and the corresponding corona inception voltages.Weather conditions have no significant impact on corona distortion.The approach proposed in [19] relies upon the observation that, for voltages substantially higher than the corona inception level, the time delay as a function of travel distance becomes linear and the steepness of the overvoltage is independent of the voltage value.The following relationship has been proposed in standards [20][21]: where S o is the original steepness of the overvoltage, S is the new steepness after the waveform travels for a distance d, and A is a constant.The constant A is a function of the line geometry only and is dependent as well on the surge polarity.Typical values are given in [20] and [21].Equation ( 4) can be used to estimate the variation with travel length of the steepness of lightning overvoltages that impact a substation [1].6.Several models have been proposed to represent insulator strings in lightning studies.The most accurate representation relies on the application of the leader progression model [1], [22][23], which can be used to account for non-standard lightning voltages.According to this model, the flashover 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 flashover occurs when the leader crosses the gap between the cross-arm and the conductor, so the total time to flashover is the sum of the corona inception time, the streamer propagation time and the leader propagation time.For details about the calculation of these times see [4], [17].Figure 2. Linear corona model.
Other models are based on the so-called integration methods and voltage-time curves.Integration methods are based on the following assumptions: (i) there is a minimum voltage V 0 that must be exceeded before any breakdown can start or continue; (ii) the subsequent time-tobreakdown is a function of both the magnitude and the time duration of the applied voltage above the minimum voltage V 0 ; (iii) there exists a unique set of constants associated with breakdown for each insulator configuration.In the most general formulation, different weights can be given to the effects of voltage magnitude and time.Although easy to use, these methods can be applied to specific geometries and voltage shapes only.Voltage-time (or time-lag) curves give the dependency of the peak voltage of the specific impulse shape on the time-to-breakdown, see Figure 3 [17].These curves are experimentally determined for a specific gap or insulator string and may be represented with empirical equations, applicable only within the range of parameters covered experimentally [1]: In practice, measurements can be affected by several factors: impulse front shape, front times of the applied standard lightning impulse, gap distance and gap geometry, polarity, internal impedance of the impulse generator (due to the predischarge currents in the gap).
7. Footing impedance modeling can be based on a nonlinear frequency-dependent circuit [24], [25].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 [17].The value of this resistance is approximated by the following expression [1], [20], [23]: 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 sufficient soil ionization.I g is given by where r is the soil resistivity (ohm-m) and E 0 is the soil ionization gradient (400 kV/m, [25]).8.A lightning stroke is represented as a current source whose parameters, as well as its polarity and multiplicity, are randomly determined according to the distribution density functions recommended in the literature [1], [7].See previous section.

IMPLEMENTATION OF THE MONTE CARLO PROCEDURE USING PARALLEL COMPUTING
This section details the most important aspects of the procedure implemented in MATLAB for calculating the LFOR of overhead transmission lines using parallel computing.The main steps are basically those implemented in previous works [3], [6].A summary of the new procedure follows.
1.The values of random parameters are estimated.
In this work, this step includes the calculation of the parameters of the lightning stroke (peak current, rise time, tail time, location of the vertical leader channel) and the phase conductor voltages.2. The last step of a return stroke is determined by using the electrogeometric model, see Figure 3. Volt-time characteristic.
Figure 4 [1].The values of the striking distances are calculated according to the following expressions: (9) where α and γ are constants that depend on the object, β is a constant that depend on the electrogeometric model and I p is the stroke peak current [1], [2].3. The previous step decides the final target of each lightning stroke: shield wire (tower or midspan point), phase conductors, ground.4. Strokes to ground are discarded, while those to the line (either to a tower or a midpsan point) are edited for including the random parameters into the ATP file of the test line model.That is, once the point of impact of each lightning stroke is known the ATP file that will simulate each stroke to the line is edited.5.The MATLAB application distributes the ATP files between the cores of the installation.6. Overvoltage calculations are performed.The only difference between models for backflash and shielding failure simulations is the node to which the current source that represents the stroke must be connected.The result of each simulation (i.e.flashover, no flashover) is reported.In case of flashover the simulation is stopped.7. The MATLAB application post processes the result of each simulation to obtain the LFOR and any other distribution function of interest (e.g. the distribution of lightning stroke peak currents that caused flashover).
The convergence of the Monte Carlo method can be checked as the simulations progress or prior to any simulation by comparing the calculated probability density functions of all random parameters to the theoretical functions.In this work, the second option is applied, so a high enough number of cases is selected to guarantee the convergence of the method prior to simulate any case.in which the current source that represents the stroke and the random parameters (i.e.lightning stroke parameters, phase conductor voltages) are embedded.Table 1 lists the random values that can be generated by the application implemented for this work [3], [6].In the present work insulator string and footing resistance parameters are not random; instead they are fully specified in the ATP overhead line template prior to the calculation of any random value.

ILLUSTRATIVE EXAMPLE Test Line
Figure 6 shows the tower design for the line tested in this paper.It is an actual 400 kV line located in Northern Algeria.As shown in the figure the line has two conductors per phase and one shield wire, whose characteristics are provided in Table 2.

Line and Lightning Stroke Parameters
A model of the test line was created using ATP capabilities and following the guidelines summarized above.
• The line is represented by nine 390-m spans plus one 10-km section as line termination at each side of the point of impact.Each span is divided into thirteen 30-m sections (for corona effect modelling) represented by means of frequency-dependent distributed-parameter models.• The towers are represented by means of the so-called multistory model presented in [26].
See also [16] and [17] for more details.• Only single-stroke negative-polarity flashes are considered.A return stroke is repre-sented by a concave waveform (see Figure 1), with parameter n = 3.24 in equations (2a) and (3).• Footing resistances are represented as nonlinear resistances modelled according to equation (7) with R o = 20 W. The value of the soil resistivity r is 100 ohm-m.• Insulator strings are represented by means of a simple model (implemented in ATP MODELS language) that duplicate a time-lag curve (see Figure 3) The length of insulator strings is 5.6 m.Assuming the critical flashover voltage (CFO) of insulator strings for negative polarity strokes is estimated as [20] CFO − = 700 ⋅ d s (10) where d S is the striking distance of insulator strings, the corresponding CFO -is 3920 kV [3].
The curve has been fitted to obtain breakdown at 2 μs with 1.58•CFO and breakdown at 3 μs with 1.36•CFO [17].The resulting coefficients for equation ( 4) are A = 2273.6E+3,B = 5448.8E+3,and m = -0.5.• Corona effect is incorporated into the line model by means of a simplified version of the circuit depicted in Figure 2, in which only the values of C 1 , R 1 , and V i are to be specified.As mentioned above, each line span is divided into sections of 30 m, and the parameters for each section (i.e.V i , R 1 , C 1 ,) are obtained according to the following equations: where r is the conductor radius, in cm, and h is the conductor height, in cm.where r is the phase equivalent radius, in cm, and h is the conductor height, in cm.K c in equation ( 11c) is adjusted to obtain a propagation model as close to equation (5) as possible.An adequate value of K c should be between 0.5 and 1.5; in this work K c = 0.6.
The following probability density functions are used to obtain random values: • Lightning stroke: Peak current magnitude, rise time and tail time are determined by assuming a log-normal distribution for all of them.Table 3 shows the values used for each parameter.All parameters are assumed independently distributed.• Stroke location: It is randomly estimated by assuming a vertical path and a uniform ground distribution of the leader.Vertical channels are uniformly distributed in a surface with a single line span length (i.e.390 m) and within a 500 m distance at each line side (see Figure 7).The electrogeometric model is represented in equation ( 9) with α = 8, γ = 0.65, and β = 1 [1].
• Phase conductor voltages: The reference angle is estimated by using a uniform distribution between 0 and 360 degrees.

Simulation Results
The LFOR of the test line was estimated by generating up to 100000 combinations of random numbers (i.e.those needed to characterize the lightning stroke and phase conductor voltages).Figure 7 shows the distribution of vertical channels within that surface after estimating the location of a few hundreds of vertical channels.
Once the coordinates of the channels were known, the electrogeometric model was applied to obtain the strokes that would impact the line and discard those to ground.Actually, less than 9000 strokes did finally impact the line (see Table 4).Figure 8 provides the distribution of lightning strokes to the line, but distinguishing between those that impacted the shield wire from those that impacted a phase conductor.
A summary of simulation results is presented in Figure 9 and Table 4.The figure provides the results obtained with two different values models of the current source that represents the lightning stroke.Although the percentage of strokes to phase conductors was not negligible, the SFFOR was in all case studies zero.From the application of the electrogeometric model, the maximum peak current magnitude of strokes to phase conductors was always below 20 kA, and none stroke could cause flashover.The percentage of strokes to the shield wire was much higher and the flashovers caused by these strokes began with peak current values above 60 kA.Consequently, the distribution of stroke currents that caused flashover was that of strokes to the shield wire.
Given the surface within which the randomly generated lightning vertical channels were distributed and assuming that the flash ground density was N g = 1 flash per km 2 and year, the generation of 100000 vertical channels corresponds to an analysis of the test line during 256410 years.Consider the case with a resistance of 400 W without including corona effect; since only 1806 strokes could cause flashovers, then the estimated LFOR for this test line is 1.806 flashovers per 100 km and year.Remember that there is a proportionality with the flash ground density, so if this density was higher than 1 the estimated LFOR of the test line should be increased in the same ratio.
As expected, the number of flashovers increases with the value of the source resistance, and reaches its maximum value for an ideal source.One can observe that the addition of corona to the line model reduces the number of flashovers, and this effect increases with the parallel resistance of the source current that represents the stroke; with an ideal source current the final value of LFOR is about 6.5% lower with corona effect, see Table 4.
The study was carried out with a multicore installation in which the ATP simulations were distributed between 50 cores.The simulation time required to carry out the entire study, including the time needed by the MATLAB application to pre-edit ATP files and post-process ATP results, and that needed by the ATP to calculate overvoltages (and decide whether the line will flashover or not) was less than 12 minutes in all case studies, see Table 4.This is a significant advance with respect to previous experiences.The work presented in [3] and [6], which can be considered a first version of the Monte Carlo approach presented here, required more than 7 hours of single CPU time to obtain the results provided in Figure 9.

DISCUSSION
• The study has been carried out with some simplifications: several parts of the transmission line model have not been represented by the most accurate models (e.g.insulator strings, corona effect), voltages induced by electric and magnetic fields of lightning channels to shield wires and phase conductors are neglected, a vertical direction is assumed for the stroke leader when it approaches earth, only single-stroke negative-polarity flashes were assumed and using an incidence model based on the electrogeometric model, which is a very simple approach for representing a very complex physical phenomenon as lightning.• Corona effect has been represented using a very simple model, which is not accurate enough for multiconductor line studies.The approach applied here was used to obtain a line model in which corona could introduce a steepness variation close to that recommended in equation ( 3); the goal is to have a line model that could be used for estimating the incoming surge [1].• The BFOR is sensitive to the coefficient of correlation between the peak current magnitude and the rise time, although the SFFOR is not; consequently, the LFOR decreases as the coefficient of correlation increases; see [3] and [10].In other words, accurate knowledge of this coefficient is an essential issue for flashover rate calculations.• Sensitivity studies can be very useful to analyze the influence of line and stroke parameters, and to determine what range of values might be of concern [3].Although the number of parameters involved in lightning calculations is very high, only some of them can be accurately specified from the line geometry.For results presented in previous works [3], it is obvious that the flashover rate increases with the peak current magnitude and decreases with the rise time, while the influence of the footing resistance is not critical for low values of the soil resistivity and low values of R o .

CONCLUSIONS
A statistical analysis is the natural approach in calculations when some system or model parameters are not well known or they are known with some uncertainties.This is a common situation in overvoltage calculations because both stress and strength are usually characterized by probability density functions from which a failure rate can be derived.When the calculations have to be obtained from the application of a time-domain software tool, such as an EMTP-like tool, a very long simulation time can be required to perform a statistical study based on the Monte Carlo method.
Although very powerful hardware platforms and very accurate software tools are currently available, simulation times of an order of hours will be usually needed.A solution for this drawback is the usage of a multicore environment.Parallel computing is a straightforward solution for parametric and statistical studies that require very long simulation times since in both cases tasks can be distributed between several cores to significantly reduce CPU time.
This paper has proposed a MATLAB application to be used in combination with the ATP package to achieve this goal.MATLAB has been selected for this Buehren for allocating ATP tasks in a multicore environment [27].
The main goal of the paper was to prove the advantages that a multicore installation can have on statistical studies carried with an EMTP-like tool.The transmission line model used for this study could have been implemented with more accurate models for some line components (i.e, insulator strings); therefore, the results derived from this study should be used with care.Not only the incidence model but other aspects needed to obtain the LFOR of an overhead transmission line have been simplified.See [3], [10] and [28] for other aspects that could have been considered in this work.

Figure 5 Figure 5 .
Figure 5 shows a diagram of the procedure implemented in MATLAB.In this work, ATP capabilities are used to edit the line model template

Figure 8 .
Figure 8. Distribution of lightning strokes to the line.

Table 3 .
Statistical Parameters of Return Strokes.

Table 4 .
Lightning Performance of the Test Line (100000 runs, 50 cores).application because none EMTP-like tools has been already developed for taking advantage of parallel computing; instead MATLAB allows users to create open environments that take advantage of high-level computing capabilities and the possibility of using parallel computing.The procedure implemented in this work uses the modules developed by M.