A three-dimensional direct lightning strike model for lightning protection of the substation

Lightning strikes to a substation happen frequently, which has a serious impact on the safety of the power system. Therefore, it is important to evaluate the effect of the direct lightning to the substation for the sake of protection. However, since the propagation and attachment of the lightning leader is a typical random process, it is hard to sufﬁciently represent the effect of the direct lightning by some deterministic models. In order to take into consideration the randomness of the lightning propagation process, in this paper, a three dimensional lightning fractal model that reproduces the tortuosity and branching of natural lightning is established. The proposed model is then applied to the effect analysis of the direct lightning to a typical 220 kV substation. Moreover, the protection effect of different distribution and height schemes of lightning rods is investigated using the proposed model. This work helps to understand the lightning strike risk that the substations face and also provides a promising way to improve


INTRODUCTION
The substations are regarded as one of the most essential components in the power system. Since the substations are normally located at the remote and open plain, the lightning may be easy to strike the substation and then cause a serious impact on the power system [1,2]. Therefore, it is important to evaluate the effect of the direct lightning to the substation for the sake of protection. However, because the substation configuration is rather complex and the real experiment is impractical, how to make a quantitative evaluation of the lightning strike risk as well as the protection performance will be a key challenge. Many researchers have developed methods to evaluate the direct lightning risk and the protection performance. Among them, some methods mainly make an engineering evaluation of the protection zone as a result of a lightning rod based on the pure theory [3][4][5], whereas some other methods mainly analyse the impact of a certain lightning current and voltage on a certain struck object [6][7][8]. However, since the propagation and attachment of the lightning leader is a typical random process, it is hard to sufficiently represent the effect of the direct lightning by such deterministic methods. It becomes a key point to make the random process of the propagation and attachment of lightning leader computable, and to consider the influence of the form of the struck object.
For that purpose, in 1980s, Niemeyer et al. proposed a simulation algorithm to reproduce the fractal growth approach of dielectric break-down channel, which was called NPW model [9]. Later, Dellera and Garbagnati studied the physical process of lightning propagation toward the earth and proposed the leader progression model (LPM) [10,11]. They analysed the initiation and propagation process of downward leaders and upward leaders of the natural lightning and proposed a simulation method to reproduce it. Since then, many lightning models were developed based on the NPW and LPM methods. Among these models, the so-called lightning fractal model that considers the tortuosity and branching nature of lightning channel progression is widely used in the simulation study of lightning strike probability [12][13][14][15][16][17]. Kawasaki and Matsuura adopted a mathematical concept of fractal dimension as a verification index for the lightning model and applied it to the simulation of a strike to a single tall building [12]. Zhang et al. adopted a 2D lightning model to calculate the strike distribution on single lightning rods of various heights [13]. Petrov et al. developed a 3D lightning fractal model to make predictions of the probability of lightning strikes to square buildings [14]. It was also applied to study the shielding failure of transmission lines [15,16]. All of these studies are based on the lightning fractal model and have been applied successfully to the evaluation of the lightning strike risk or protection performance.
However, the existing lightning models cannot be directly used in the evaluation of the lightning protection effect of the substation. Specifically, the existing models only involve struck objects with a simple geometry such as a long transmission line or a square building. They choose several inception positions of the lightning upward leader at the tip points of the simple structure, but it is much more difficult to choose inception points in a substation that has so many complex structures. Moreover, the lightning fractal model has not been used to evaluate the immunity performance of a complex structure group with different voltage levels and lightning protection requirements, such as a substation. The traditional model gives a lightning strike probability analysis of different attachments, but does not judge whether the equipment is damaged by the direct lightning.
In this paper, a 3D lightning strike model for complex structures will be proposed. Based on the proposed model, the lightning strike probability and direct lightning current for different equipment can be obtained. Moreover, the proposed model will be applied to the evaluation of the direct lightning strike risk and the protection effect for a 220 kV substation. Compared with the previous work, the proposed model is mainly improved in three aspects: 1. The lightning model in a 3D space is established, in which the inception points of the upward leader are selected quantitatively, rather than qualitatively, by calculating the electric field on the envelope surface of earthed structures. 2. The box fractal dimension of the 3D lightning is calculated and used as a verification index for quantifying and reproducing the tortuosity and branch characteristic of natural lightning flashes. 3. Combining the basic insulation level of different facilities in the substation with their lightning strike distance, a simplified criterion for evaluation of the direct lightning damage is adopted. The lightning strike probability and lightning current of different equipment are obtained based on the proposed model.
The remainder of this paper is organized as follows. Section 2 describes the modelling method of the lightning strike process and the direct lightning risk assessment for a group of complex structures by the model. In Section 3, the sensitivity of model parameters is discussed and the proposed model is applied to the evaluation of the direct lightning risk and the protection performance for a 220 kV substation. Finally, Section 4 includes the summary and conclusion for this paper.

Lightning strike model
In the LPM, the downward leader initiates from thunderclouds and causes the upward leader to initiate from a structure tip or corner [10]. They attract each other in the space until the lightning final-jump occurs. Actually, the progression of lightning leaders takes place in steps [25]. So, the lightning strike process can be calculated by iteration as shown in Figure 1. To compute the iteration process, space is divided by small cubes with a side of step length L, and the lattice points are defined as "developed points" inside the lightning leaders and "candidate points" near the lightning channels, respectively. Table 1 shows the parameters of the lightning strike model, and the sensitivity of model parameters will be discussed in next chapter. It should be noted that E T is not the mean cloud-to-ground electric field equal to U T /H T , but the electric field near ground [19]. The critical upward leader inception electric field of a transmission line is calculated by [24] E cl = 3000 k(1 + 0.03∕ where k is the surface roughness factor, δ is the relative air density, r is the radius of the transmission line. At each iteration step, the direction of lightning leader progression depends on the spatial electric field. The probability of the candidate point turning into a new developed point could be expressed as [21]  L The minimum elongation of the leader at each iteration step 1 m [ 20] where E b is the critical background electric field, E m→n is the average electric field between the mth developed point and the nth candidate point, N d is the number of developed points, N c is the number of candidate points, and η is an uncertain positive number called development probability exponent. The average electric field E m→n along the side of the small cube can be calculated by where U m is the space potential of the mth developed point, U n is the space potential of the nth candidate point. As the leader develops along the diagonal of the cube, the denominator of Equation (3) changes accordingly.
The electrical potential U m in the downward leader can be calculated by where L m is the length of the lightning leader from the thundercloud to the mth developed point. The space potential including U n at each point outside the lightning leaders is calculated with the finite difference method of Laplace's equation [17] The successive over-relaxation (SOR) iterative algorithm is used to solve Laplace's equation under the boundary condition U T , U s and E T . The space potential at the (i+1)th iteration step can be calculated by [26] where U x,y,z is the space potential at the point (x, y, z), ω is the SOR factor. The factor ω is selected as 1.95 in the subsequent simulation.
To quantify and reproduce the tortuosity and branch characteristic of natural lightning flashes, fractal dimension is used as a verification index of simulated lightning [12]. As a mathematical concept, fractal dimension has many different expressions. The most widely used fractal dimension, box dimension, can be defined as [15] where D is the box fractal dimension, a is the side length of measurement box, N(a) is the number of measurement boxes to cover the whole measured object. The development probability exponent η in Equation (2) has a big influence on the shape of lightning branches as shown in Figure 2. Therefore, in this paper, the value of η is adjusted to control the fractal dimension of simulated lightning. According to the observation results of the natural lightning, the fractal dimension of lightning is usually between 1.1 and 1.4 [12,16,21]. As shown in Figure 3, η = 5 is suitable for reproducing natural lightning where the fractal dimension is about 1.2.
A concept of "hot-spots" is used to solve the problem that thousands of points on the structure surface may satisfy the electric field condition E cs or E cl simultaneously [15]. The hotspots are preset inception points of upward leaders with a steep electric field, and usually located at a structure tip with large curvature. For each hot-spot, the lightning strike procedure can be described as shown in

Boundary potential calculation
The proposed model can be used to simulate the path of lightning channels in the huge cloud-to-ground space. However, for a less demand on computer memories and CPU resources, there is no need to consider the lightning fractal branches until the distance of the downward leader head from the earthed structure is less than the maximum strike distance of the structure [14]. Thus, the simulation space can be set in a much smaller space. As shown in Figure 5, the simulation space has a side boundary with the potential gradient E T . Its upper boundary potential is decided by the equipotential thundercloud plane and a vertical lightning without fractal from thunderclouds to the upper boundary. The fractal downward leader of the proposed model starts from the intersection of the vertical leader and the upper boundary. The boundary condition of the simulation space is fixed for a simplified calculation, which

2.3
Direct lightning risk assessment

Lightning strike probability
The lightning strike probability describes the chance of one facility in all structures being attached by downward lightning leader, which can also be used to describe the ability of a lightning rod to attract lightning. In this paper, to obtain the lightning strike probability for a facility among all the facilities of the complex structures group, the Monte Carlo idea is adopted. In the method, the lightning strike probability for a facility can be represented by the relative strike frequency, which can be expressed as where R is the relative lightning strike probability among all facilities, P L is the real lightning strike probability of each facility, F is the lightning strike frequency of each facility in simulation tests.
Equation (8) holds only if simulation tests are done so many times. Therefore, hundreds of simulation tests for each case will be performed to get the lightning strike probability.

2.3.2
Simplified direct lightning damage criterion This paper mainly focuses on the construction of the3D lightning strike model that can reflect the propagation and attachment process of the lightning strike. After the lightning strike occurs, a simplified criterion is adopted for evaluation of the direct lightning damage. It is an empirical and incomplete criterion, which does not consider some factors such as the lightning current reflection and flashovers.
When the lightning current flows through the attachment into the ground, the struck object is subjected to a lightning voltage. The criterion assumes that the struck object will not be damaged unless the lightning current or voltage on the equipment is much larger than its basic insulation level (BIL) when shield failure occurs [27]. The relationship between the minimum lightning current and the BIL of equipment is given by [28] where I is the lightning strike current that devices can withstand, U BIL is the BIL of equipment, Z char is the surge impedance of conductor.
The strike distance refers to the distance between the head of the downward leader and the attachment when the final-jump occurs. The electrical geometric model requires that lightning strike distance is interrelated to lightning strike current The critical safe strike distance corresponding to different BIL magnitude, which could be expressed as [29] S = 10I 0.65 s (10) where S is the lightning strike distance, I s is the lightning strike current magnitude.
According to Equations (9) and (10), it can be judged whether the direct lightning damages its attachment in a simulation test. Taking Z char as 300 Ω [30], the critical safe strike distance with relation to BIL of a vertical wire remote from earth is shown in Figure 6. To be noted, this equation is only a rough calculation for assessing the possibility of damage in this paper. Actually, the direct lightning risk requires a more complicated assessment process.

RESULTS AND DISCUSSION
In this section, the lightning strike model will be applied to the direct lightning protection of the substation. A model of the typical 220 kV substation in China is shown in Figure 7, and the side and top views of the substation with its dimensions are given in Figure 8. For simplicity, this substation just includes the most necessary and important electrical equipment as possible lightning strike receptors. The BIL of substation equipment is set according to the design manual of State Grid.
Here are some necessary instructions about the 220 kV substation model: 1. The location and size of the substation apparatus are finetuned according to the real 220 kV substation, so that the tip is located at a 1 m grid point. 2. Towers, equipment shells and edges of the building with lightning strap installed are regarded as ground. 3. Due to the lightning's high velocity in the range of 10 5 to 10 6 m⋅s −1 [20], the power frequency voltage is not considered. The three phases of the line are all set to the maximum phase voltage as a simplified condition.

Simulation results
In reality, the lightning may come from all directions in the sky. But the number of starting points of the downward leader has to be limited for the realization of simulation. There are three As shown in Figure 9, the electric field of the substation envelope surface under thunderclouds is calculated with the finite difference method before simulation. Twenty points that have the relatively large and steep electric field are selected as the hotspots.
The distribution and height of lightning rods in the substation are designed by the rolling sphere method [31]. The two different configurations of the distribution of lightning rods, the fourrod configuration and the five-rod configuration, are shown in Figure 10. The height of all rods is 30 m and will be changed in subsequent simulations.
The lightning strike probability should be obtained from a large number of simulation tests according to the Monte Carlo idea. In some past studies, this number has ranged from 100 to 1000 [13][14][15]. In this paper, 300 simulation tests are performed for each case. These simulations are carried out on a high performance computing platform using 380 CPU cores of 2.6 GHz. The average time for a single simulation is 2.3 h.
According to different lightning attachments, the simulation results can be divided into three cases: 1) lightning strike to the rod (LTR), as shown in Figure 11(a); 2) lightning strike to the ground or the side boundary (LTG), as shown in Figure 11(b) and; 3) lightning strike to the equipment or structure (LTE), as shown in Figure 11(c,d). The lightning strike probability of these three cases under different lightning rods distributions is shown in Table 2. In addition, some model parameters are changed and their sensitivity is discussed in the following subsection.

Sensitivity analysis for the model parameters
The values of the lightning strike model parameters in Table 1 are mainly derived from others' experiments or observations. However, the results of experiments and observations are different. For example, the leader inception field E cs may be reduced to 400 kV/m and the internal leader field E d may be reduced to 6 kV/m [32,33]. In this subsection, the lightning strike probability and mean strike distance of 220 kV line, transformer and  towers are used as the observation objects to explore the sensitivity of the obtained results to some of the model parameters.
The parameter E cs varies from 400 to 600 kV/m and the corresponding results are shown in Figure 12(a,b). It can be seen that the lightning strike probability and strike distance of towers are highly relevant to E cs . As the leader inception field increases, the lightning strike probability and strike distance of towers are both reduced, The breakdown field E g seems to have little effect on the results as shown in Figure 13(a,b). It may indicate that the randomness of lightning strike probability and strike distance is mainly due to the random propagation process of the lightning leader before the final-jump.
The internal leader field E d varies from 5 to 15 kV/m and the corresponding results are shown in Figure 14(a,b). As the internal field increases, the potential of the downward leader tip decreases. It is more difficult for the initiation of the upward leader, and it is easier for lightning to strike the boundary and ground. Thus the strike probability and strike distance decreases.
As shown in Figure 15(a,b), the boundary field E T caused by thunderclouds has little effect on the results. It should be noted that the influence of the thundercloud potential U T on the boundary potential is similar to that of E d and E T , therefore U T will not be discussed in this paper.

Vulnerability analysis for the equipment
For lightning protection, it is important to determine which equipment of the substation is vulnerable to direct lightning. By drawing the Pareto chart of lightning strike frequency of different equipment under the no-rod scheme as shown in Figure 16, we find that the 220 kV line, the transformer and the towers have a much bigger chance to be struck than other equipment or structures. But the damage rate of the 220 kV line and the transformer is very different. By a histogram distribution analysis as shown in Figure 17, more than a quarter of lightning strikes to the 220 kV line exceed its current limit, whereas about 90% of lightning strikes to the transformer are below its BIL. The reason may be that the electric field distortion on the surface of the wire with positive charges is larger than the grounded equipment, which is easier for initiation of the upward leader.

Analysis of protection effect of different lightning rods distributions
The number of lightning rods and their distribution are the main factors in lightning protection if the rod height does not change. However, it is difficult to compare the protection effect of different lightning rod distribution schemes by theoretical calculation. In this paper, the direct lightning protection effect of the four-rod configuration and the five-rod configuration introduced above is quantitatively evaluated by the proposed method. The lightning strike probability of vulnerable equipment under different distributions of lightning rods is shown in Figure 18. It can be seen that the five-rod configuration reduces the probability of lightning strike to vulnerable equipment much greater than the four-rod configuration. Only about 10% of lightning strikes are not incepted by the five lightning rods.

FIGURE 16
The counts of lightning strike to different attachments without lightning rods protection. The line means the Pareto cumulative frequency of them However, the mean strike distance of vulnerable equipment under the four-rod scheme and the five-rod configuration are approximately the same, which is about 10 m less than that under the no-rod configuration, as shown in Figure 19. It shows that an additional lightning rod cannot further reduce the lightning current of equipment significantly.  The detailed counts distribution of lightning strike to rods is shown in Figure 20. Compared with the strike distance distribution of the four-rod configuration, the five-rod configuration reduces the lower limit of strike distance from 30 to 23 m, which means a larger incepted current range and a better protection effect for equipment. It also explains the high probability of the five rods incepting the lightning in Table 2.

Influence of the rod height on strike probability
In this subsection, the rod height of the four-rod configuration is adjusted and the change of lightning strike probability of the lightning rod and other attachments is obtained. As shown in Figure 21, the probability of strike to lightning rods The lightning strike probability of the lightning rod and equipment under different lightning rod heights tends to increase rapidly first and then stabilize gradually or even decrease with the increase of the height of lightning rods. Accordingly, the probability of lightning strike to equipment shows a generally opposite trend.
The trend of strike probability of lightning rods below 45 m is consistent with the theory of EGM. However, the strike probability of the rod height of 50 m is abnormally reduced. To confirm whether the taller rod is inefficient in the smaller simulation space of 120 × 120 × 120 m, a series of simulation tests are carried out in a 240 × 240 × 240 m space as shown in Figure 22. The y coordinate of the starting point is 100, 120 and 140 m. As shown in Figure 23, the strike probability of lightning rods of 45 m and above has been improved, which is more consistent with theoretical expectations. The results show that there is a threshold of the height of the struck object for effective calculation of the strike probability. The height threshold is related to the size of the simulation space. When the height of the rod is under the threshold, the lightning strike probability increases with the height of the rod and this trend is basically consistent in a larger simulation space. When the rod height exceeds the threshold, it is necessary to expand the simulation space to obtain effective results.

CONCLUSION
In this paper, a 3D lightning strike model based on LPM is proposed and applied to the assessment of the direct lightning protection performance for a 220 kV substation in China. By calculating the electric field on the envelope surface of the substation under thunderclouds and combining the BIL of equipment with the strike distance, the lightning strike probability and damage rate for different equipment are obtained. What is more, the sensitivity of the obtained results to model parameters is explored.
The results show that the leader inception field and internal leader field are parameters with higher sensitivity. Under the selected model parameters, the 220 kV power line is the most vulnerable equipment when lightning strike occurs. When four lightning rods are distributed in the corner of the substation, a rod added in the centre can greatly reduce the risk of lightning strikes to vulnerable equipment. What is more, the results show that there is a threshold of the height of the struck object for effective calculation of the strike probability, which is related to the size of the simulation space. The lightning strike probability increases with the rod height under the height threshold and this trend is basically consistent in a larger simulation space.