RESEARCH ARTICLE Fractal characteristics of simulated electrical discharges

A two dimensional stochastic dielectric breakdown model was utilized to study the fractal properties of simulated lightning discharges and surface discharges. The fractal dimension of the simulated growth patterns varied depending on the cell configuration chosen for the breakdown. The inclusion of the cells in the diagonals with reference to the growth site produced less branched trees with smaller fractal dimensions. The production of branches in the electrical tree growth depends highly on ‘ η ', which is the exponent of the breakdown probability distribution. For small η values, highly complex tree patterns with many branches were observed. By controlling the value of η , growth patterns similar to the experimental observations could be produced. The average fractal dimension estimated through the Correlation Function method when η =1 for simulated lightning discharges and surface discharges were 1.56±0.01 and 1.68±0.01 respectively. When η >4 the growth patterns effectively lose their fractal structure and became a curve with dimension 1. Keywords: Electrical treeing, lightning, simulation, surface discharges doi:10.4038/jnsfsr.v36i2.145 Journal of the National Science Foundation of Sri Lanka 36 (2) 137-143


INTRODUCTION
Stochastic models for dielectric breakdown processes provide a possible theoretical foundation to describe the tortuous structure of lightning discharges in the atmosphere as well as the development of surface discharges on insulating materials. In 1984, a simple two dimensional stochastic model was produced to support the experimental observations of branching gas discharges 1 . The model was developed as a growth pattern, which finds the new growth sites according to the potential in the nearest neighbours by solving the discrete Laplace equation. The basic assumption is that the growth probability of conducting structures depends on the local electric field. They showed that their model naturally leads to fractal structures with dimension 1.75 ± 0.02 for simulated discharge patterns. Others Previous studies 3 analyzed a set of lightning photographs and concluded that the average fractal dimension of lightning discharges is 1.34 ± 0.05, and built a simple dielectric breakdown model in two dimensions to explain the quantitative fractal behaviour observed in lightning photographic images. The two dimension stochastic model has been extended to simulate lightning discharges in three dimensions and calculated the average fractal dimension in 3D as well as the average fractal dimension in the 2D projections 4 . They reported a value ≈1.51 as the fractal dimension for simulated 3D lightning patterns and ≈1.34 for the vertical projection of the same. Their observations agree with previously published results only at higher value of the exponent η in the breakdown probability distribution.

Fractal characteristics of simulated electrical discharges
In a recent study, another group presented an improved leader progression model to simulate lightning channels 5 . They calculated the path of lightning channels and their progression under the condition of finite-charge distribution in the thundercloud and developed numerical simulation by computing the charge distribution in a square lattice. They also defined the charge distribution using the Poisson's equation, and for each time step, finite-difference method was applied to solve the Poisson's equation numerically. Their estimation on fractal dimension varied in the range from 1.1 to 1.4.
Most theoretical models that deal with the propagation of lightning generated electromagnetic fields do not take the fractal nature of lightning into account. A recent study 6 , which has investigated the tortuosity and the branching of lightning channels, has shown that the fractal nature of lightning can significantly modify the radiated waveforms and increase the high frequency content. Incorporating methods that could simulate fractal characteristics similar to the experimental measurements make the lightning propagation models more realistic.
Dielectric materials are often used as passive components in electrical appliances. The condition of the dielectric material determines the reliability and lifetime of electrical systems. It has been reported that the fractal dimension of electrical breakdowns decreases with the decrease in the dielectric constant of the insulator 7 . Therefore, it is important to investigate the influence of the local electric field which is a proxy for the dielectric constant on the fractal development.
The main objective of this work was to study the influence of initial cell configuration on the fractal geometry of simulated lightning discharges and surface discharges, and estimate the value of η in the breakdown probability distribution that produces fractal structures similar to the experimental observations. The work presented here is based on the Laplacian growth model which is similar to the original model developed in previous studies 1 . Several fractal dimension estimation algorithms were used in interpreting the simulated growth patterns, and the results were compared with the published values for natural dielectric breakdown processes.

METHODS AND MATERIALS
Dielectric Breakdown model: The simulation was developed to generate growth patterns of lightning discharges and surface discharges in two dimensions. The initial charge configurations for both processes are shown in Figure 1. These two configurations act as the initial boundary conditions for the subsequent development of the growth patterns. The black dots represent the cells having positive potential ( =1) and the grey dots represent the cells having negative potential ( =0).
Both, lightning and surface discharge configurations were simulated on a 250 x 250 square lattice. For lightning configuration, a negative charge was placed above the layer of positive charges at a distance of 250 lattice units. In surface discharge configuration, the negative charge was placed at the centre of the lattice and a circle of positive charges were placed around the centre charge at a radial distance of 250 lattice units. The model works as follows. The pattern grows stepwise from negative potential ( =0) to positive potential ( =1). The electric potential ' ' in each cell was calculated by solving the Laplace equation numerically over the 2D grid (Equation 1) subjected to the initial boundary conditions ( Figure 1). In the finite-difference approach, for a square grid, central difference approximation of the second derivative yields, (2) where (x, y) is an arbitrary cell coordinate and h is the distance between adjacent lattice points. Thus, by using potential of the nearest neighbours in x-direction and y-direction, potential at an arbitrary cell (x, y) can be calculated; i.e. the potential at any grid point is the mean potential of its nearest neighbours. To increase the speed of convergence, the system of equations generated through Equation 2 was solved by using successive overrelaxation method 8 as shown in Equation 3.
Here, the index n and n+1 refer to the previous and current iteration cycles of the system respectively. The over-relaxation parameter λ, has the effect of amplifying each step towards the final answer when λ>0. In this work, a value of 0.5 was used for λ. When the residuals of all nodes were below a predefined tolerance, the iterative process was terminated which was found to be approximately the size of the lattice.
Step by step, the algorithm searches the possible lattice points which can be attached to the growing pattern. At each step, the Laplace equation was solved numerically and one lattice point was added to the pattern as a new bond with zero potential. Each time a new point is added, it alters the potential of the unoccupied cells in the vicinity. The new point is chosen from all possible adjacent cells subjected to the probability distribution P i given in equation 4. where the index i refers to the list of adjacent cells and n is the total number of cells. Basically, P i acts as a weight function to choose cells according to their potential. After each step, the potential of all unoccupied cells was recalculated to introduce the effect due to the newly added cell. This process was repeated until the pattern reached the positive potential (cell boundary). The effect of the exponent η on the growth patterns are discussed later.
It should be noted that the discharge development process discussed here, represents the movement of the stepped leaders for lightning discharges and streamers for surface discharges. If one wishes to develop the growth model further, for example to introduce dart leaders to follow the same general path created by the stepped leaders in the lightning configuration, or, distortions of the electric field due to space charges distributed along the streamer channels in the surface discharge configuration, then after completing the initial growth cycle, the Laplace equation must be replaced with a Poisson equation for the subsequent growth cycles. Since the first growth cycle deposits charges along the leader channels or the streamer channels, subsequent growth cycles will be automatically attracted to the old path, thereby reducing the fractal nature.
Estimation of fractal dimension: Three fractal dimension estimating methods were used to evaluate the dimension of individual growth patterns. The first method is the popular Box Counting method 9 . In this method, the number of boxes 'n(r)' inside which, at least one sample point lies were counted for different box sizes (r). The Box Counting dimension D1 is calculated using equation 5.
The second method is a variant of the Box Counting method called 'Sandbox' method 2 . In the Sandbox method a square box of size 'L' is formed on the pattern and the mass of the tree found within the box is evaluated. The mass can be evaluated by counting the number of lattice points within the box. Average mass M(L) is obtained for different box sizes 'L'. Fractal dimension D2 is the exponent that expresses the scaling of the mass with its size. The third method, called Correlation Function method, gives a statistical value for fractal dimension based on pair-wise distance between points. Basically, Correlation dimension can be calculated between each pair of points in a set of 'N' number of points. In this method, mass is equal to the average of the number of pairs where pair distance is less than the defined length 'r'. The mass can be calculated using correlation integral 10 , where (X i -X j ) is the pair-wise distance and Θ is given by, The Correlation dimension D3 can be calculated mathematically by using, 3 ( ) D C r r ∝ ......... (8) Generally, in all three methods, fractal dimension can be found by calculating the gradient of a double log (log-log) plot through the least-square fitting procedure.

Accuracy of fractal dimension algorithms
Three fractal dimension estimating techniques were utilized to extract the fractal dimension of the simulated discharge patterns. The implemented algorithms were tested by applying on well known shapes before extracting the fractal dimension of simulated discharge patterns. Two Euclidean shapes and two fractal shapes, where the theoretical fractal dimensions are known, were used as known shapes. Figures 2a and 2b show the two fractal shapes used in this work, known as Sierpinski Gasket and Sierpinski Carpet, respectively. Table 1 shows the estimated results for each of the three fractal estimating methods against their theoretical values. When the fitting errors were examined by minimizing the χ 2 , it was noted that the statistical errors in the determination of the fractal dimension was in the order of ±0.001. According to the results shown in Table 1, all fractal dimension estimation methods have produced the fractal dimension within ±0.07 of the theoretical values. Although the Box Counting algorithm has performed quite well here, it should be noted that the algorithm has a number of practical limitations in estimating dimensions of complex patterns, especially at high embedding dimensions 9 .

Effect due to cell configuration
The speed of electrical tree growth and appearance of the branches in the growing pattern depend on the selected cell configuration, which defines the possible growth sites at each step. There are three possible cell configurations that can be selected to simulate lightning and surface discharges, namely, cross configuration, diagonal configuration and space filling configuration. For these three cell configurations, growth patterns were simulated by keeping η=1. Fractal dimension of the simulated patterns were evaluated using the Box Counting, Sandbox and Correlation Function methods. Figures 3a, 3b and 3c show three typical simulated growth patterns for lightning discharges. The insets show the cell configuration used. Fractal dimension of simulated surface discharges were also found by using the same three fractal dimension estimating methods. Figures 3d, 3e and 3f show the three typical simulated growth patterns for surface discharges with the insets showing the cell configurations used.
In Table 2, estimated mean fractal dimensions and their errors for three cell configurations for both lightning and surface discharges are shown. For each configuration, 25 growth patterns were generated randomly to estimate the mean and error. The error in the mean was evaluated as σ/√n where σ is the standard deviation and n is the number of growth patterns. It can be seen that the selection of the diagonal points from the nearest neighbours lead to simulated growth patterns with less spatially separated branches contributing to a smaller fractal dimension (see cross configuration and diagonal configuration). In general, the box Counting method yields a significantly smaller fractal dimension compared to the Sandbox and Correlation Function methods. This is due to the fact that the Box Counting technique emphasized the geometrical structure of the fractal pattern ignoring the underlying measure.    4a, 4b and 4c show the dependency on η for simulated lightning discharges for cell configuration 1 for η=0.5, 1.0 and 1.5 respectively. When η increases, the 'dense' nature of the tree reduces rapidly and less branching was observed. Figures 4d, 4e and 4f show the same for surface discharge configurations for η=1.0, 2.0 and 4.0 respectively. In the surface discharge configuration, when η increases, the symmetrical development of the pattern was completely abandoned and the growth pattern tended to produce a branch or two aligned in random orientation growing from the centre of the lattice to the boundary. i.e., the growth pattern preferred to stay on course with the initial direction which was chosen randomly at the set off. This implies that, in actual situations with insulators, if there are surface defects which favour the initial discharge development, the electrical breakdown will always occur in the same path if the dielectric constant of the insulator is high enough to reduce the fractal growth. This behaviour can seriously affect the insulator and degrade the performance and increase ageing. Dulan Amarasinghe & Upul Sonnadara  The variation of fractal dimension with the exponent η for lightning and surface discharges are shown in Figure 5 for the Correlation Function method. When η>4, the growth patterns effectively lose their fractal structure and become a curve with dimension 1. For lightning configuration, when η<1, the growth pattern is highly dense, deviating from any experimentally observed discharge patterns. It should be noted that when η=0, as expected, the fractal dimension moves towards 2 which is the dimension of a plane. As shown in the figure, the relationship between the fractal dimension and η can be expressed by a 3 rd degree polynomial. Although the same growth model was used, the simulated lightning patterns lose their fractal structure faster than the simulated surface discharges. This is due to the influence of the initial charge configuration on development of the tree patterns.
Only few published results are available on fractal dimension of natural breakdown trajectories for lightning and surface discharges. The available data suggest that the fractal dimension for natural lightning and surface discharges are in the order of ≈1.34 and ≈1.70 respectively. Thus, for cell configuration 1, the dielectric breakdown model discussed in this study can reproduce discharge patterns with a similar fractal structure when η ∼ 1.6 and η ∼ 0.9 for lightning and surface discharges respectively ( Figure 6). That is, surface discharges tend to produce highly complex, tortuous paths compared to natural lightning discharges. The tortuous paths of lightning discharges are responsible for enhancing the high frequency end of the lightning generated EM spectrum. Less tortuous surface discharge channels are responsible for the breakdown of solid dielectrics. Thus, the effect of tortuosity on the lightning generated EM fields and the breakdown of solid dielectrics can be studied using the model presented here. When η =1 simulated lightning and surface discharges produced fractal dimension values of 1.56 ± 0.01 and 1.68 ± 0.01 respectively.

Execution time and growth
One of the main problems in the simulation is the execution time which increases exponentially with the increase in the lattice dimension. In the present study, this was avoided by limiting the subsequent iterations to a grid of 40x40 lattice points centred on each newly added growth site. However, this may lead to patterns growing asymmetrically especially when simulating surface discharges in large lattice structures. The execution time depends not only on the lattice dimension but also on the over-relaxation parameter 'λ' given in Equation 3. The possible 'λ' lies between 1 and 0 and it controls the speed of convergence. It was found that a value of approximately 0.5 seems to work well.
It was noted that the points added on to the existing pattern at later stages favour the development of branches. This behaviour was seen both for lightning discharges and surface discharges. Thus, a cell already embedded in the pattern has a lesser probability of branching at the later stages due to screening from the existing cells.
In the interpretation of the results presented here it must the kept in mind that the stochastic model discussed above depends on the 2D Laplace fields to determine the electric potential on cells of a 2D lattice (without considering the effects due to the dimension perpendicular to the plane of the lattice). Past studies have shown that 3D Laplace fields tend to produce fractal patterns with slightly higher dimensions compared to 2D Laplace fields for the development of surface discharges 11 . For lightning discharges, studies carried out with the 3D stochastic models have shown that, compared to 2D fractal dimension, 2D projections of 3D fractal patterns show a 10-13% reduction in fractal dimension 4 . Although the 2D dielectric breakdown model discussed in this study can be easily extended to 3D, the heavy computational time required to develop fractal structures in 3D introduces a severe constraint on the development of growth patterns even for reasonably sized lattice structures such as 250×250.

CONCLUSION
A 2D stochastic model based on Laplacian growth was used successfully to simulate lightning discharges as well as surface discharges. A strong dependency on the fractal nature of the growing patterns was seen on the chosen cell configurations. The selection of diagonal points from the nearest neighbours simulates growth patterns of lightning as well as surface discharges with lesser branches contributing to a smaller fractal dimension.
The most popular method to compute fractal dimension of a given structure is the Correlation Function method, which estimates dimension based on the pair wise distances between points on the pattern. However, since the magnitude of the fractal dimension may depend on the fractal dimension estimation method, in addition to the Correlation Function method, Sandbox method and another popular method, Box Counting method, both of which can be implemented easily, were used in this study. The differences of fractal dimension values reported in these methods are related to the emphasis on point wise and bulk size dimension. Thus, one should be cautious when comparing fractal dimension estimates based on different techniques.
Fractal dimensions of simulated discharge patterns in all cell configurations show strong dependence on the exponent η. By changing the η value, simulated discharge patterns, where the fractal dimension is closer to the experimental values, can be developed. That is, the appearance of the branches in the simulated tree patterns is controlled by the exponent. When η increases the 'dense' nature of the trees are reduced rapidly. That is, patterns effectively lose their fractal structure and become a curve with dimension 1 when η>4 .
It has been shown that the model can be used to generate growth patterns with similar fractal structures to natural lightning and surface discharges when η∼ 1.6 and η ∼ 0.9 respectively. Further work is necessary and is in progress to compare the growth patterns generated through the dielectric breakdown model with models that can generate similar fractal structures using diffusion limited aggregation.