Investigation of Arc Trajectory Characteristics of Composite Insulators Based on Multifield Coupling Analysis

This paper proposes an arc propagation model based on electric and thermal field coupling to investigate the arc propagation mechanism around composite insulators. The charge simulation method reduces the computational complexity of the electric field calculation and simulates the leakage current density and charge motion during the arc propagation. The finite difference time domain method describes the continuous variation in the thermal field generated by the leakage current density and the arc close to the insulator. Different types of heat transfer are analyzed to evaluate the arc injection and dissipation energies. The stochastic arc trajectory property is modeled by the Markov process based on electric and thermal field coupling instead of the conventional deterministic flashover model. The modified model simulates three typical arc trajectories from ignition to extinction to analyze the mutual effects of arc trajectory and arc energy. The results show that the arc trajectories that travel away from the insulator and along the leakage distance are more likely to be extinguished during propagation than the arc trajectory bridging most sheds. The consecutively rising leakage current and relatively low energy dissipation during the arc propagation process are two dominant factors that cause flashover. The model explains the mechanism of the phenomena in which the probability of arc jumping between sheds increases when the insulator surface pollution level becomes severe, which leads to a decrease in the average flashover voltage.


I. INTRODUCTION
Insulators are used to provide electrical insulation and mechanical support for transmission lines [1]. Arcs are likely to ignite close to high-voltage (HV) electrodes and propagate in the air close to the insulator surface [2]- [4]. Arcs initiate when the instantaneous electric field exceeds the threshold value and ionizes air [5]- [7]. Consecutive arcs are inclined to formalize a conductive path from an insulator HV electrode to a ground electrode and lead to flashover.
Many researchers have studied arc propagation models on the insulator surface. Obenaus and Neumarker introduced the classic static arc propagation and flashover model in 1958 [8]. Sundararajan et al.  The associate editor coordinating the review of this manuscript and approving it for publication was Yanzheng Zhu . gation and improved the model with dynamic parameters [9]. Ibrahim applied this dynamic model to washed HV insulators and reduced the flashover probability during the insulator washing process [10]. Both static and dynamic arc propagation models assume that the arc trajectory follows the leakage distance of the insulator. However, laboratory experiments indicate that arc trajectories around the insulator could be different even when the supply voltage, surface contamination level, and environmental factor remain the same. This paper provides a stochastic model to describe the process of arc from ignition to extinction by modeling charge motion under the influence of electric and thermal fields in the air.
The charge simulation method (CSM) is utilized to calculate the electric field distribution and the value of charges around an arc simultaneously. H. Singer et al. were the first to apply CSM to calculate the electric field for HV equipment in 1974 [11], [12]. Takuma improved the CSM model by considering charge density on the interface between different materials to model the capacitive-resistive field [13], [14]. El-Kishky compared the electric field distribution along the leakage distance for porcelain and composite insulators [15]- [17]. Most existing research applied the CSM to calculate the steady-state field without the rearrangement of the charge positions, while this paper utilizes the CSM to analyze the continuously redistributed electric field and space charges during arc propagation. The probability of arc growth directions depends on the instantaneous value of the electric field, which ionizes the air and moves charges to formalize the arc trajectory.
The thermal field around an arc is calculated by the finite difference time domain (FDTD) method to simulate heat conduction, convection, and radiation processes during arc propagation. Doufene presented the heat dissipation model on the polluted insulator surface to illustrate that the insulator surface temperature increases with increasing electric field strength [18]. Tatematsu used FDTD-based simulation to investigate the breakdown characteristics of a long air gap [19], [20]. Existing research demonstrates that the FDTD method is suitable to analyze the field variation process numerically in the time domain. However, the previous analysis of the arc heat transfer process did not comprehensively explain the influence of the relative fluid velocity of arc. This paper simulates the arc heat conduction and convection processes simultaneously by the FDTD based on the characteristics of the arc propagation.
The Markov process is used to describe the stochastic characteristics of arc propagation. Markov proposed the theoretical model in 1906. Many researchers have modified the Markov model by introducing the hidden chain to predict the behavior of invisible variables based on the visible results [21]. This paper models the arc propagation steps as states in the Markov process and calculate the transition probability matrix based on the present electric field distribution.
The purpose of this paper is to propose an arc propagation model based on multifield coupling to investigate the stochastic property of arc trajectories. The model presents arc energy injection and dissipation at each step of propagation to elaborate the mechanism of arc ignition and extinction.

II. MODEL SCHEMATIC AND FIELD CALCULATION
The schematic of the FXBW4-220/100 composite insulator model is shown in Figure 1. The supply voltage is set to 220 kV at the HV electrode of the composite insulator. The environment temperature is 298 K, the air pressure is the standard atmospheric pressure (101.325 kPa), the air humidity is 74% according to local weather conditions, and the wind speed is 0 m/s. Equivalent salt deposit density (ESDD) is used to describe the pollution degree on the insulator surface. The ratio T /B of the pollutant on the shed top surface to that on the shed bottom surface is defined as ESDD Top /ESDD Bottom . The flashover voltage increases with decreasing pollution degree ratio, and the range of the ratio is from 0.1 to 1 [22]. Therefore, the pollution degree ratio is selected as 1 to simulate severe contamination conditions with low flashover voltage. The ESDD value is 0.1 mg/cm 2 under the influence of temperature and humidity. Water particles in the air are not considered in the model, since the humidity is less than the air saturated humidity [23].

A. ELECTRIC FIELD CALCULATION MODEL
The Poisson equation describes the quasi-static electric field distribution. The field satisfies the Dirichlet and Neumann boundary conditions.
where ϕ is the potential distribution in the field domain, ρ c is the charge density, and ε 0 is the permittivity of the insulating material. CSM is applied to numerically calculate the instantaneous electric field around the arc. Ring charges are utilized to simulate the axis-symmetric geometry of the composite insulator and reduce the computational complexity [24]- [26] (Figure 2). Line charges are utilized to model the arc trajectory. Point charges represent the free charges in air and describe the electron avalanche at the arc leader. The contour points on the interfaces of the HV electrode are shown in Figure 2.   The boundary condition on the interface between the electrode and air is calculated in (2): The boundary condition on the interface between the electrode and insulator is calculated in (3): where Q E represents ring charges in the HV electrode; Q I and Q A are ring charges in the insulator and air near the insulator surface, respectively; and ϕ E is the electrical potential of the HV electrode. m I and m A is the number of charges in the insulator and air, m E is the number of charges in the electrode. The arrangements of simulated charges and contour points on the interface between the air and insulator are shown in Figure 3.
To maintain potential continuity on the interface between the insulating material and air (ϕ I = ϕ A ), the Dirichlet boundary condition is expressed as To maintain flux density continuity on the interface between the insulating material and air (ε I E I − ε A E A = σ ), the Neumann boundary condition is expressed as Charge density σ is calculated by the leakage current density on the insulator surface: where ε I and ε A are the permittivity of the insulator material and air, respectively. S(i) are small surface areas surrounding contour points. R(i) is the resistance of the insulator surface between two contour points. C(i) is the capacitance of the insulator surface between two contour points. ω is the radical frequency of the AC power system. Contour points and charges on the arc trajectory are shown in Figure 4.
The boundary condition on the arc trajectory is expressed as where Q line represents line charges in the arc trajectory, and m arc represents the number of line charges in the arc. 6132 VOLUME 8, 2020

B. THERMAL FIELD CALCULATION MODEL
Charge motion and arc propagation satisfy the principle of fluid flow. Thermal transfer between the air and arc consists of heat conduction, convection, and radiation.
Heat conduction is the form of heat transfer on the interface between the arc and air with direct contact. The heat conduction partial differential equation (PDE) and boundary conditions are expressed as where T is the temperature in the field domain, and t is time. ρ, c and λ are the fluid density, specific heat capacity, and thermal conductivity of the mediums, respectively. v is the velocity field, and is the internal heat generated by the arc and leakage current on the insulator surface.
To solve the PDE, the thermal field close to the arc is discretized by cubic voxels with the side length of h( Figure 5). Each vertex of the cubic voxel represents a discrete point.
The PDE (8) is discretized with FDTD method as where u, v, and w are the X-, Y-, and Z -components of the fluid velocity, respectively. Heat convection is the heat transfer process with relative motion between the arc and air. An arc with an internal heat source flows through the air and insulator surface.
The convection motion between the arc and air is described by mass and momentum conservation equations expressed as where ρ is the fluid density, and p is the atmospheric pressure. u, v, and w are the X-, Y-, and Z -components of the fluid velocity, respectively. η is the viscous stress tensor. The thermal field distribution is calculated based on (8) and (11). The heat transfer coefficient α i on the interface between the insulator material and the arc is calculated as The heat transfer coefficient α a on the interface between the air and arc is calculated as The heat transfer energy caused by conduction and convection is calculated as where W heat is the heat transfer energy. T (t) is the temperature of the arc. T a and T i are the temperature of the air and insulating material respectively. λ a and λ i are the thermal conductivity of the air and insulating material respectively. S(l, t) is the heat transfer area as a function of time and arc length. l a and l i represent the lengths of the arc in air and on the insulator surface, respectively. t 0 is the period from heat transfer initiation to the moment of calculation.
Heat radiation plays a significant role in the heat transfer of an arc. Heat radiation is the process of arc internal energy decrease caused by emitting radiant energy. The heat ray follows the laws of reflection and refraction. Thermal radiation energy W radiation is calculated as where ξ is the ratio of the thermal radiation subject to the thermal radiation of the black body. T is the internal thermal temperature of the arc, and δ is the Stefan-Boltzmann constant. t 0 is the period from radiation initiation to the moment of calculation.

III. STOCHASTIC ARC PROPAGATION MODEL A. ARC DIRECTION BASED ON ELECTRIC FIELD
The stochastic arc propagation model applies the Markov process to investigate the random property of the arc. The Markov process describes a sequence of possible events in which the probability of each event depends on the state attained in the previous event.
Arc Position n = Arc Position n−1 P ij = (Arc Position n−1 P ij )P ij = Arc Position n−2 P 2 ij = (Arc Position n−2 P ij )P 2 ij = Arc Position 0 P n ij (16) where n is the number of states. P is the transition probability matrix, which contains the probabilities that convert to the next state. i represents the current state, and j points to the next possible states. The arc previous position, current position, and next position represent three arbitrary states in the sequence of the Markov process. The state transfer function is shown in Figure 6.
In Figure 6, the arc trajectory moves from the previous position to the current position with transition probability P 12 and moves from the current position to the next position with probability P 23 . The transition probabilities are calculated by the electric field distribution around the arc trajectory in (17). A random number is generated to determine the exact direction for the next step of arc propagation. Therefore, the arc trajectory could be different even when the electric field distribution remains the same, which is consistent with the stochastic property of arc propagation. However, the distribution pattern of arc trajectories can be found when arc propagation processes are repeated multiple times. The arc ignition criterion is that the maximum electric field exceeds the threshold value of air breakdown.
where E is the electric field strength summation of all possible positions with E > E c , and E c (2.1 kV/mm) is the RMS value of the threshold field. θ is the step function [25]. The flashover probability is calculated by where N A represents the number of arcs that complete the flashover, and N T is the total number of arc propagation processes [27].

B. ARC ENERGY BASED ON THERMAL FIELD
The arc energy consists of the energies that the arc absorbs and dissipates during propagation. The arc injection energy is calculated as where u arc (x, t) and i leakage (x, t) represent the potential and the leakage current on the arc trajectory at an arbitrary time during propagation, respectively. t 0 is the period from radiation initiation to the moment of calculation. l is the length of the arc trajectory. Arc energy dissipation is the sum of heat conduction, convection, and radiation. (20) The arc energy is calculated as

C. STOCHASTIC ARC PROPAGATION MODEL BASED ON THE FIELD COUPLING EFFECT
The arc propagation model includes the arc direction determined by the electric field and the arc energy based on the thermal field. The electric and thermal fields are coupled by the instantaneous field calculation and heat transfer simulation. The coupling effect of electric and thermal fields is shown in the flow chart ( Figure 8). From Figure 8, it is observed that the instantaneous electric and thermal fields are calculated at each state in the Markov process. The electric field affects the thermal field by altering the arc propagation direction and the shape of the arc trajectory. The thermal field influences the electric field by determining the arc next state based on the arc energy. The coupling effects of the electric and thermal fields are combined with the Markov process to analyze the stochastic property of arc propagation.

IV. RESULTS AND DISCUSSIONS A. ELECTRIC AND THERMAL FIELD DISTRIBUTIONS
Instantaneous electric and thermal fields are calculated in this section to create the arc trajectory and obtain the arc energy during propagation. Arc energy is utilized to explain arc extinction and flashover criteria. The arc trajectory and  thermal field distribution are shown in Figure 9 when the arc trajectory reaches Shed 21.
Electric field distributions along the insulator leakage distance are shown in Figure 10 when the arc approaches different sheds. Figure 10 shows the electric field distribution along the arc trajectory and the remaining leakage distance of the insulator during arc propagation. The maximum electric field is located at the arc leader. The arc trajectory formalizes a conductive VOLUME 8, 2020 path, and it shortens the remaining leakage distance of the composite insulator. The potential is subjected to reduced leakage distance and leads to an increase of the electric field. The flashover probability increases significantly when the arc trajectory covers 3/4 of the insulator leakage distance because the continuously increasing electric field vectors to all the directions around the arc leader exceed the dielectric strength of air.

B. ARC ENERGY OF DIFFERENT TRAJECTORIES
Arc bridging large sheds is considered the phenomena of arc jumping between sheds. Arc trajectory 1 is defined as the arc trajectory on which the occurrences of arc bridging to adjacent sheds are less than 15% during the propagation, representing the arcs that travel along the insulator surface. Arc trajectory 2 is defined as the arc trajectory on which the occurrences of arc bridging to adjacent sheds are more than 15% but there is no occurrence of arc jumping over four consecutive sheds, representing the arcs that jump between sheds. Arc trajectory 3 is defined as the arc trajectory with the occurrence of arc jumping over four consecutive sheds, representing the arcs that travel away from the insulator. These three types of arc trajectories include all arc propagation phenomena. Figure 11 shows three typical examples of each arc trajectory category. However, the arc trajectory could be slightly different even when they belong to the same category due to the stochastic property of the arc.
Sheds 1-4 are selected in the first instance to analyze the arc properties during the propagation process.   It is observed in Figure 12 that Trajectories 1 and 2 achieve higher injection energy than Trajectory 3 because Trajectories 1 and 2 reach the surface of Shed 4. The leakage currents of Trajectories 1 and 2 increase significantly due to the decrease in the remaining surface resistance without the arc. Trajectory 3 does not reach any insulator surface after Shed 1, so the leakage current remains the same during the further arc propagation process. Figure 13 presents the case in which the dissipation energy of Trajectory 1 is higher than those of the other two trajectories because Trajectory 1 travels the longest distance in the air and has the largest heat transfer area.
Arc energy comparisons of three arc trajectories are shown in Figure 14. The arc energy of Trajectory 3 decreases to zero. Therefore, Trajectory 3 stops propagation and is extinguished. Figures 15 and 16 show the arc injection and dissipation energy of Trajectories 1 and 2, respectively, when they approach Shed 18.   Figure 15 shows that the injection energies of arc Trajectories 1 and 2 reach the same value due to the same remaining contamination layer resistance at Shed 18. Figure 16 illustrates that arc dissipation energy of Trajectory 1 increases significantly faster than that of Trajectory 2, because Trajectory 1 jumps more insulator sheds than Trajectory 2, so the heat transfer area and propagation time of Trajectory 1 are greater than those of Trajectory 2.
The arc energies of Trajectories 1 and 2 are shown in Figure 17 when they approach Shed 18. Figure 17 shows that the arc energy of Trajectory 1 decreases to zero. Trajectory 1 stops propagation and is extinguished, while Trajectory 2 continues to propagate and causes flashover.

C. ARC ENERGY WITH DIFFERENT ESDD VALUES
Arc Trajectory 1 is selected to explain the mechanism of the fact that flashover probability increases with increasing insulator surface pollution degree. The injection energies with different ESDD values are shown in Figure 18.    Figure 18 indicates that the increase in arc injection energy with large ESDD values is higher due to the significant growth in the leakage current. Therefore, it is concluded that an arc trajectory with high injection energy has a high VOLUME 8, 2020 probability of finishing the conductive path from the HV electrode to the ground electrode before extinction.
The insulator flashover probabilities of three arc trajectories are shown in Figure 19. The arc propagation processes are repeated 437 times at each ESDD value to obtain the adequate number of three typical arc trajectories and fit the curve. Figure 19 shows that the arc trajectory 1 and 3 have the probability to cause flashover due to the stochastic property of the arc, the arc trajectories could be slightly different even they belong to the same category.
The total flashover probability is divided into three parts, P flashover = P surface P 1 + P jump P 2 + P away P 3 (22) where P surface represents the proportion of arcs that travel along the insulator surface out of the total number of arcs. P jump represents the proportion of arcs that jump between insulator sheds out of the total number of arcs. P away represents the proportion of arcs that travel along the insulator surface out of the total number of arcs. P 1 , P 2 , and P 3 are the flashover probabilities of the three arc trajectories, respectively. The ESDD value has significant positive effects on the increase of P 2 and P jump . The increase in P 2 is caused by the increase in the leakage current and arc injection energy. The increase in P jump results from the increasing proportion of the vertical electric field vector to the ground insulator during the arc propagation process.
The average flashover voltage has been calculated to compare with the experiment results of insulator Type D in [28], which has a similar dimension to the insulator model used in this paper. Figure 20 demonstrates that the simulation model has similar results to the experiments. The error between simulation and experiment flashover voltages is 7.17%. The error is caused by the stochastic property of arc propagation. The flashover experiments were conducted three times to obtain the flashover voltage, and the average difference between maximum and minimum flashover voltage under different ESDD values is 6.55% [28], while the simulations were repeated for 437 times to calculate the average flashover voltage. It is difficult to repeat the experiments for the adequate number of times to eliminate the effects of the arc stochastic property.

V. CONCLUSIONS
This paper proposes a stochastic arc propagation model to describe the mechanism and behavior of an arc around a composite insulator. The combined effects of the electric and thermal fields on the arc trajectories are analyzed by the Markov model and arc energy variation during the arc propagation process. The conclusions of the paper are as follows: 1. The electric and thermal fields are coupled by the instantaneous field calculation and arc energy model. The electric field determines the directions of arc propagation and the shape of the arc trajectory, which is the heat source of the thermal field. The thermal field influences the instantaneous electric field by computing the arc energy and determining the next state of the arc as continuing propagation or extinction. 2. The coupling effects of electric and thermal fields are combined with the Markov model to analyze the stochastic property of arc propagation instead of the traditional flashover model. The arc trajectories are categorized into three types in order to investigate the effects of arc trajectory patterns on the flashover. The arc trajectory jumping between sheds has a flashover probability of 96.3% when ESDD value is 1 mg/cm 2 because the leakage current increases consecutively, and the energy dissipation is relatively low during the propagation process. 3. With the analysis of the Markov transition matrix and arc trajectories distribution, it is concluded that an arc trajectory with a large ESDD value has a higher probability of bridging sheds than an arc with a small ESDD value because the vertical electric field vector to the ground insulator proportionally increases with increasing ESDD values.