Development of a New Bi-Arc Dynamic Numerical Model for Modeling AC Flashover Processes of EHV Ice-Covered Insulators

This paper presents the development of a new bi-arc dynamic numerical model for predicting AC critical flashover voltage (FOV) of ice-covered extra-high voltage (EHV) insulators. The proposed model is based on a generic calculation algorithm coupled with commercial finite element method software designed to solve the Obenaus/Rizk model. The proposed model allows one to implement the Nottingham and Mayr approaches and compare the results obtained as a function of the arcing distance, the freezing water conductivity, and the initial arc length. The validation of the model demonstrated high accuracy in predicting the FOV of ice-covered post-type insulators and its capability to simulate the interaction of the two partial arcs during the flashover process. In particular, the results showed that the Nottingham approach is sensibly more accurate than the Mayr one, especially in simulating the dynamic behavior of the partial arcs during the flashover process. Based on the encouraging results obtained, a multi-arc calculation algorithm was proposed using the bi-arc dynamic numerical model as a basis. The basic idea, which consists in dividing the multi-arc model in several bi-arc modules, was not implemented and validated but will serve as a promising concept for future work.


Introduction
Flashover of line and post insulators due to atmospheric ice accretion still constitutes a cause of failure of overhead electrical power transmission systems in cold climate regions.From several reports from several utilities, the flashover phenomenon of high voltage (HV) insulators can occur during ice accretion, or after, under warming conditions [1][2][3][4][5][6][7][8].Several studies performed on this topic demonstrate that the flashover process is an extremely complex phenomenon, which depends on the interaction between the partial arcs established along the air gaps, the condition of the ice surface, the environmental conditions, and the arcing distance of the insulator [1][2][3][4][5][6][7][8][9][10][11].
In the last several decades, the advancement of knowledge flashover process on ab ice-surface has led to the development of several single arc mathematical models, both static [12][13][14] and dynamic [15][16][17][18], dedicated to predicting flashover voltage (FOV) of short ice-covered insulators.These mathematical models were based on the Obenaus/Rizk approach developed for polluted insulators.This approach is an equivalent electrical circuit where a single arc is in series with a residual electrical resistance [19,20].In the case of polluted insulators, the residual resistance is calculated from the Wilkins analytical expression, which takes into account the arc root in contact with the ice surface [21].However, the Wilkins formulation, initially developed for single arc model, limits the applicability of the static and dynamic predictive models to insulators having an arcing distance lower than 1 m [12][13][14][15][16][17][18].
To avoid such a limitation, some authors have proposed an improvement of the original Wilkins formulation in order to determine the residual resistance when several partial arcs and, consequently, several arc roots are in contact with the ice surface [14,18].This improved formulation has been implemented in a mathematical static multi-arc model and successfully applied to full-scale extra-high voltage (EHV) post insulators [14].More recently, this new formulation has been implemented in a dynamic mathematical multi-arc model using the same approach and is described in more detail in the next section [18].Both multi-arc models can predict the critical flashover voltage, but the dynamic model can also predict leakage current and arc velocity during the flashover process.However, the use of Wilkins formulation requires a well-defined uniform ice layer that can only be obtained under severe icing conditions [14].Under these specific conditions, the insulator is totally bridged by the ice deposit that can then be modelled as a half cylinder [12][13][14][15][16][17][18].Consequently, the calculation of the residual resistance using the analytical formulation remains the main limitation of the current static and dynamic mathematical models used to predict the FOV of ice-covered insulators.However, such residual resistance formulation is more difficult to use in the case of non-uniform conductive layers with complex geometries.
To deal with such geometry-related problems, new numerical predictive models using E-field calculation tools have been proposed in recent years by the authors [22][23][24].The initial numerical model proposed was based on the Obenaus/Rizk single arc model, which was solved using the finite element method (FEM) to compute the E-field distribution between the arc root and the ground electrode, the leakage current as well as the residual resistance of the ice surface for a defined applied voltage.The accuracy of this new predictive model was significantly improved for FOV results obtained with single arc mathematical models without previously associated geometrical limitations thanks to FEM.These results lead to the development of a numerical static bi-arc model extending the numerical single arc model to longer insulators of arcing distance up to 2.02 m [22].The accuracy of the results obtained with this bi-arc model also showed improvement compared to mathematical models.However, these static numerical models do not take into account the implementation of the arc velocity displacement and do not provide information about the temporal evolution of the arc at the ice surface as well as the leakage current.Such limitations of the numerical static models were solved by the recent development of a generic numerical single arc model in which the Bondiou-Gallimberti arc velocity criterion was implemented to handle the dynamic behavior of the single arc [23,24].This model was validated with both experimental results obtained from ice-covered and polluted insulators with an improvement in the accuracy of FOV compared to mathematical and numerical static models.
Based on the promising results obtained with the numerical dynamic single arc model developed in recent work, the authors decided to extend this model to the original numerical static bi-arc model presented in detail in [22].For that, two different approaches proposed by Nottingham [25] and Mayr [26] to simulate the partial arc were implemented in the bi-arc numerical dynamic model.The results were then compared and validated in terms of FOV, leakage current, and evolving partial arc velocities.
The proposed bi-arc dynamic numerical models proposed in this paper represent a great improvement over actual mathematical multi-arc models, which are confined to simple geometries without providing accurate and convivial numerical tools for outdoor insulator dimensioning.

The Obenaus/Rizk Single Arc Flashover Model
The first mathematical model of flashover was proposed by Obenaus and dedicated to polluted insulators [19].This simple model, illustrated in Figure 1, is the basis of most of the flashover models developed in the last four decades for both the polluted and ice-covered insulators.For Obenaus, flashover can be modelled as a single arc of length x in series with a conductive layer of length L−x, where L is the total arcing distance of the insulator.The conductive layer of residual resistance R(x) is used to model the humidified polluted layer at the surface of the insulator or the water film formed on the insulator ice surface during the melting period.The corresponding circuit equation for the Obenaus model can then be expressed as follows [19]: where Vm is the applied voltage (V); x is the local arc length (m); Im is the leakage current (A); R(x) is the residual resistance of ice layer (Ω); Ve is the electrode voltage drop.The latter can be neglected under AC, where Earc is the electric field along the arc, as initially proposed by Nottingham and expressed as follows [25]: with A and n being the arc constants.Equations ( 1) and ( 2) describe the flashover process for DC applied voltage.For AC, however, Rizk observed that the arc extinguishes and re-ignites twice in each cycle [15].The arc re-ignition occurs when the AC applied voltage Vm reaches a specific value, the arc re-ignition condition, which can be expressed as follows [20]: where k and b are the arc re-ignition constants, and Vm and Im are respectively the peak value of the AC applied voltage (V) and the corresponding peak value of current (A).
In order to solve Equation (1), the residual resistance R(x) must be determined.For mathematical flashover models, two approaches have been principally used: the open or AR model proposed by Rumeli and the mathematical formulation proposed by Wilkins.Rumeli proposed that the surface geometry of a real insulator is equivalent to a surface in 2D, called the AR or open model, whose length is equal to the leakage path of the given real insulator, as explained in detail in [27].The advantage of this model is to provide a simple representation of the complex geometry of the uniformly polluted layer covering the insulator.However, the calculation of the pollution resistance of this open model requires that the density current distribution is uniform.In this context, the presence of the arc root at the surface of the conductive layer cannot be taken into account, resulting in a significant discrepancy in the calculation of R(x), as demonstrated in a previous study [27].On the other hand, Wilkins proposed a new formulation, which takes into account the effect of the density current line constriction at the arc root position.Wilkins considered the arc root to be circular so that the arc is cylindrical [21].In this way, using the theory of conjugate functions, the calculation of the polluted surface resistance is reduced to a two-dimensional (2D) Laplacian field problem.Taking into account the presence of the arc root influences the residual resistance R(x) value as demonstrated in [27].However, the problem with the Wilkins formulation is that it can only be applied to rectangular pollution layer geometries and is not representative of real polluted insulator shapes.In fact, this formulation has been principally applied in the modeling of ice-covered insulators that present a rectangular conductive layer in the case of a uniform ice layer [12][13][14][15][16][17][18], as illustrated by Figure 2. In most ice-covered insulator flashover models, the simplified For Obenaus, flashover can be modelled as a single arc of length x in series with a conductive layer of length L−x, where L is the total arcing distance of the insulator.The conductive layer of residual resistance R(x) is used to model the humidified polluted layer at the surface of the insulator or the water film formed on the insulator ice surface during the melting period.The corresponding circuit equation for the Obenaus model can then be expressed as follows [19]: where V m is the applied voltage (V); x is the local arc length (m); I m is the leakage current (A); R(x) is the residual resistance of ice layer (Ω); V e is the electrode voltage drop.The latter can be neglected under AC, where E arc is the electric field along the arc, as initially proposed by Nottingham and expressed as follows [25]: with A and n being the arc constants.Equations ( 1) and ( 2) describe the flashover process for DC applied voltage.For AC, however, Rizk observed that the arc extinguishes and re-ignites twice in each cycle [15].The arc re-ignition occurs when the AC applied voltage V m reaches a specific value, the arc re-ignition condition, which can be expressed as follows [20]: where k and b are the arc re-ignition constants, and V m and I m are respectively the peak value of the AC applied voltage (V) and the corresponding peak value of current (A).
In order to solve Equation (1), the residual resistance R(x) must be determined.For mathematical flashover models, two approaches have been principally used: the open or AR model proposed by Rumeli and the mathematical formulation proposed by Wilkins.Rumeli proposed that the surface geometry of a real insulator is equivalent to a surface in 2D, called the AR or open model, whose length is equal to the leakage path of the given real insulator, as explained in detail in [27].The advantage of this model is to provide a simple representation of the complex geometry of the uniformly polluted layer covering the insulator.However, the calculation of the pollution resistance of this open model requires that the density current distribution is uniform.In this context, the presence of the arc root at the surface of the conductive layer cannot be taken into account, resulting in a significant discrepancy in the calculation of R(x), as demonstrated in a previous study [27].On the other hand, Wilkins proposed a new formulation, which takes into account the effect of the density current line constriction at the arc root position.Wilkins considered the arc root to be circular so that the arc is cylindrical [21].In this way, using the theory of conjugate functions, the calculation of the polluted surface resistance is reduced to a two-dimensional (2D) Laplacian field problem.Taking into account the presence of the arc root influences the residual resistance R(x) value as demonstrated in [27].However, the problem with the Wilkins formulation is that it can only be applied to rectangular pollution layer geometries and is not representative of real polluted insulator shapes.In fact, this formulation has been principally applied in the modeling of ice-covered insulators that present a rectangular conductive layer in the case of a uniform ice layer [12][13][14][15][16][17][18], as illustrated by Figure 2. In most ice-covered insulator flashover models, the simplified Wilkins formulation given by Equation ( 4) is used with the assumption that w L, where w is the width of the ice layer, and L the total arcing distance of the insulator: where x is the arc length, γ e is the surface conductivity of the conductive layer, and r is the arc root radius.
In previous studies, the arc-root radius r and the surface conductivity γ e were established to be [12]: ) where I m is the leakage current, B a constant dependent on the nature of the applied voltage and conductive layer, σ is the conductivity of the freezing water, and α and β are constants depending on the applied voltage polarity [12][13][14].
Energies 2016, 9, xFOR PEER REVIEW 4 of 22 Wilkins formulation given by Equation ( 4) is used with the assumption that w << L, where w is the width of the ice layer, and L the total arcing distance of the insulator: where x is the arc length, γe is the surface conductivity of the conductive layer, and r is the arc root radius.
In previous studies, the arc-root radius r and the surface conductivity γe were established to be [12]: where Im is the leakage current, B a constant dependent on the nature of the applied voltage and conductive layer, σ is the conductivity of the freezing water, and α and β are constants depending on the applied voltage polarity [12][13][14].Based on the Obenaus/Rizk single arc model, several mathematical models, static or dynamic, were developed, predicting with a good accuracy the FOV of ice-covered insulators [12,13,15].However, this single arc model could only be applied to ice-covered insulators with arcing distance lower than 1 m.As reported in several experimental studies on short suspended insulator strings or small porcelain post insulators, the flashover process of small ice-covered insulators seems to be attributed to the propagation of only one electrical arc on the ice surface [4,5].For longer insulators like those used on EHV electrical networks, several electrical partial arcs on the ice can be involved in the flashover process [11,14,18].In that situation, a single arc flashover model cannot be used anymore and some modifications need to be done, as presented in the next sections.

The Principle of the Static Mathematical Multi-Arcs Model
When the arcing distance of an insulator is greater than 1 m, as in the case of EHV insulators, the dynamic ice accretion process combined with the presence of a high electric field leads to the formation of several air gaps (ice-free sections of the insulator).These air gaps combined with the presence of a water film generated by the melting of the ice surface lead to a non-uniform Based on the Obenaus/Rizk single arc model, several mathematical models, static or dynamic, were developed, predicting with a good accuracy the FOV of ice-covered insulators [12,13,15].However, this single arc model could only be applied to ice-covered insulators with arcing distance lower than 1 m.As reported in several experimental studies on short suspended insulator strings or small porcelain post insulators, the flashover process of small ice-covered insulators seems to be attributed to the propagation of only one electrical arc on the ice surface [4,5].For longer insulators like those used on EHV electrical networks, several electrical partial arcs on the ice can be involved in the flashover process [11,14,18].In that situation, a single arc flashover model cannot be used anymore and some modifications need to be done, as presented in the next sections.

The Principle of the Static Mathematical Multi-Arcs Model
When the arcing distance of an insulator is greater than 1 m, as in the case of EHV insulators, the dynamic ice accretion process combined with the presence of a high electric field leads to the formation of several air gaps (ice-free sections of the insulator).These air gaps combined with the presence of a water film generated by the melting of the ice surface lead to a non-uniform distribution of the electric field along the insulator [9,10].Under melting conditions, the increase of the electric field along the air gaps can initiate partial arcs along each air gap, as illustrated in Figure 3. Depending on the applied voltage amplitude and the ice surface condition, these partial arcs will start propagating along the ice-covered insulator, and will join together to form a complete flashover [9,10].The presence of several partial arcs will give rise to several arc roots on the ice surface.This specific condition generates a non-uniform current distribution along the melted ice surface, which directly affects the calculation of the residual resistance of the ice layer.
In order to be able to model this multi-arc flashover process in the case of ice-covered EHV insulators, some authors have proposed extending the single-arc static model, proposed by Obenaus/Rizk, to multi-arcs (Figure 4a) by dividing the ice layer covering the EHV insulator in several parts in order to create several single arc models in series (Figure 4b) [14,18].Moreover, in the case of a partial arc having two roots in contact with the ice surface (Arc 2 on Figure 4a), it must be divided into two partial arcs of equal length (Arc'2 and Arc"2 on Figure 4b) in order to create two single arc models.Hence, using this strategy, an ice-covered insulator (Figure 4a) with three partial arcs can be modelled using four single arc models in series (Figure 4b), which are governed by Equations ( 1) to (3).
Energies 2016, 9, xFOR PEER REVIEW 5 of 22 distribution of the electric field along the insulator [9,10].Under melting conditions, the increase of the electric field along the air gaps can initiate partial arcs along each air gap, as illustrated in Figure 3. Depending on the applied voltage amplitude and the ice surface condition, these partial arcs will start propagating along the ice-covered insulator, and will join together to form a complete flashover [9,10].The presence of several partial arcs will give rise to several arc roots on the ice surface.This specific condition generates a non-uniform current distribution along the melted ice surface, which directly affects the calculation of the residual resistance of the ice layer.
In order to be able to model this multi-arc flashover process in the case of ice-covered EHV insulators, some authors have proposed extending the single-arc static model, proposed by Obenaus/Rizk, to multi-arcs (Figure 4a) by dividing the ice layer covering the EHV insulator in several parts in order to create several single arc models in series (Figure 4b) [14,18].Moreover, in the case of a partial arc having two roots in contact with the ice surface (Arc 2 on Figure 4a), it must be divided into two partial arcs of equal length (Arc'2 and Arc''2 on Figure 4b) in order to create two single arc models.Hence, using this strategy, an ice-covered insulator (Figure 4a) with three partial arcs can be modelled using four single arc models in series (Figure 4b), which are governed by Equations ( 1) to ( 3).GND distribution of the electric field along the insulator [9,10].Under melting conditions, the increase of the electric field along the air gaps can initiate partial arcs along each air gap, as illustrated in Figure 3. Depending on the applied voltage amplitude and the ice surface condition, these partial arcs will start propagating along the ice-covered insulator, and will join together to form a complete flashover [9,10].The presence of several partial arcs will give rise to several arc roots on the ice surface.This specific condition generates a non-uniform current distribution along the melted ice surface, which directly affects the calculation of the residual resistance of the ice layer.
In order to be able to model this multi-arc flashover process in the case of ice-covered EHV insulators, some authors have proposed extending the single-arc static model, proposed by Obenaus/Rizk, to multi-arcs (Figure 4a) by dividing the ice layer covering the EHV insulator in several parts in order to create several single arc models in series (Figure 4b) [14,18].Moreover, in the case of a partial arc having two roots in contact with the ice surface (Arc 2 on Figure 4a), it must be divided into two partial arcs of equal length (Arc'2 and Arc''2 on Figure 4b) in order to create two single arc models.Hence, using this strategy, an ice-covered insulator (Figure 4a) with three partial arcs can be modelled using four single arc models in series (Figure 4b), which are governed by Equations ( 1) to ( 3).GND This arc splitting modelling requires modifying the expression of the residual resistance R(x) given by Equation ( 4) for the one arc root in contact with the ice surface.In a general way, if there are N arcs, N of them with one-root and N of them with two roots in contact with the ice (N = N + N ), then the total residual resistance R(x) of the multi-arc model can be expressed as follows [14]: where L is the total arcing distance of the insulator, and x is the sum of all arc lengths.Despite the relatively good accuracy obtained with this static multi-arc model, as demonstrated in [14], it does not allow the behavior of partial arcs at the ice surface to be modeled or the arc velocity and leakage current evolution during the flashover process to be computed.For that purpose, a new dynamic multi-arcs model has been proposed recently, as presented in the next section.

The Principle of Dynamic Mathematical Multi-Arcs Models
The dynamic mathematical multi-arc model proposed by [18] uses exactly the same approach as the static multi-arc model does, and is shown in Figure 4.In fact, the authors use a single arc dynamic mathematical model, presented in detail in [10], to solve the multi-arc problem by splitting of several single arc models, as illustrated in Figure 4b.The residual resistance R(x) is calculated using the same analytical formulation as used for the multi-arc static model, given by Equation (7).
The main goal here is to propose a new arc velocity criterion.The proposed arc velocity criterion is an experimental criterion based on high-speed measurements of the propagation of the upper and lower partial arcs during the flashover process.From the results obtained, the authors have determined that the average arc velocity v(x) can be expressed as follows [18]: where x(m) is the length of the partial arc.The proposed multi-arc dynamic model allows for predictions of the FOV of EHV ice-covered insulators with good accuracy (Figure 3), with an arcing distance up to 4.17 m.However, this model presents two main limitations.The first limitation comes from the experimental arc velocity criterion given by Equation ( 8), which does not permit any simulation of the velocity of each partial arc on the ice surface.Indeed, each partial arc propagates differently depending on its initial length and position along the insulator.Generally, the upper partial arc located at the top of the insulator in contact with the HV electrode is faster than the lower partial arc in contact with the ground electrode [14,18].The second limitation of this model, which is the same as that of the static multi-arcs model, is inherent to the residual resistance R(x) calculation given by Equation (7).Since this equation can only be applied to rectangular geometries as those obtained for uniform ice layers, and it cannot cover non-uniform ice layers as those obtained in the case of insulators equipped with alternated sheds or booster sheds [4,5].

The Single Arc Dynamic Numerical Model
Despite the fact that both static and dynamic mathematical models applied to multi-arcs flashover modeling provide satisfactory results in terms of FOV, they present limitations, as mentioned previously.To extend the applicability of predictive models in terms of geometry and dynamic parameter computation, the authors have developed in recent years a new dynamic single arc model to predict the FOV of both polluted and ice-covered insulators.This model, presented in detail in [23,24], uses the finite element method (FEM) commercial software Comsol Multiphysics to compute the leakage current I m and the residual resistance R(x) related to the Obenaus/Rizk model described by Equations ( 1)-(3).The numerical approach based on FEM allows for the development of a versatile model, which is independent of the geometry of the problem and, consequently, can be applied to different conductive layer geometries as those of polluted insulators using the AR model or ice-covered insulators [24].The FEM is coupled with a generic calculation algorithm implemented in Matlab (2017), which allows for the determination of the evolution of the leakage current I m and the arc velocity v(t) during the flashover process as well as the average E-field E avg in the conductive layer.In order to implement dynamic calculations, the algorithm uses the arc propagation criterion proposed by Hampton [28] given by Equation ( 9) as well as the arc velocity criterion proposed by Gallimberti [29] given by Equation (10).These equations are expressed as: where E avg is average E-field inside the pollution layer, and E arc is the electric field along the arc given by Equation ( 2), and as: where Q i (t) is the required charge to induce the arc propagation, which can be expressed as follows [29]: where C i is the capacitance between the arc root and the ground, and V ap the potential at the arc root.
In Equation (11), the capacitance C i can be calculated using a spherical approximation [29] as follows: with where L is the insulator arcing distance, x the arc length, and r the arc root radius given by Equation (5).
As mentioned previously, this single arc dynamic numerical model has been successfully applied to both polluted and ice-covered insulators presenting short arcing distance.As for longer insulators, where two and more partial arcs are involved in the flashover process, the authors decided to improve their single arc numerical model in order to take into account the presence of several partial arcs at the surface of the ice-covered insulator.

The Principle of the Bi-Arc Dynamic Numerical Model
Figure 5 illustrates the modelling principle used in the case of two partial arcs in contact with the ice surface.As observed in Figure 5, during experimental tests, each partial arc is initiated at the HV and ground electrode, respectively; consequently, each of them featuring an arc root in contact with the ice surface.Due to the presence of the two arc roots on the ice surface, the distribution of the current density is modified as well as the residual resistance R(x) between the two arc roots [22].However, this modification can be easily handled using FEM, where the influence of each arc root is taken into account without having to divide the ice layer into two parts, as done for the two mathematical multi-arc models.
The proposed bi-arc dynamic numerical model is based on the static bi-arc numerical model developed by the authors dedicated to ice-covered insulators and presented in detail in [22].To this static model, the authors have applied the dynamic conditions used in their numerical single arc model described previously to enable it to compute the evolution of the leakage current and the velocity of each partial arc during the flashover process.Hence, both the Hampton and Gallimberti criteria described by Equations ( 9) and (10) were used, respectively.In the FEM model, two arc roots in contact with the ice surface was considered as an equipotential surface of radius r given by Equation (5).The voltage boundary condition applied to each arc root was calculated as follows: Energies 2018, 11, 2792 8 of 22 where V m is the peak value of the AC applied voltage; x 1 and x 2 are the length of the Partial Arcs 1 and 2, respectively (Figure 5); I m is the leakage current; R(x) is the residual resistance of the ice layer; A and n are the AC arc constants.
Finally, as proposed in [14], the condition of arc re-ignition given by Equation ( 3) must be modified in order to take into account the presence of the two arcs one of which propagates downwards (Arc 1) and the other upwards (Arc 2).This is done by the following equation: where x 1 and x 2 are the length of Partial Arcs 1 and 2, respectively, and k 1 , k 2 , and b are AC parameters given by [12].
Energies 2016, 9, xFOR PEER REVIEW 8 of 22 in contact with the ice surface was considered as an equipotential surface of radius r given by Equation (5).The voltage boundary condition applied to each arc root was calculated as follows: where Vm is the peak value of the AC applied voltage; x1 and x2 are the length of the Partial Arcs 1 and 2, respectively (Figure 5); Im is the leakage current; R(x) is the residual resistance of the ice layer; A and n are the AC arc constants.Finally, as proposed in [14], the condition of arc re-ignition given by Equation ( 3) must be modified in order to take into account the presence of the two arcs one of which propagates downwards (Arc 1) and the other upwards (Arc 2).This is done by the following equation: where x1 and x2 are the length of Partial Arcs 1 and 2, respectively, and k1, k2, and b are AC parameters given by [12].
(a) (b) Figure 6 shows the calculation algorithm of the dynamic numerical bi-arc models.This algorithm is based on the calculation algorithm developed for the dynamic numerical single arc presented in detail in [23] with the addition of the new conditions expressed by Equations ( 14)-( 16).The problem is solved in 2D using the electric current module and the stationary solver of Comsol Multiphysics, where only the conductive surface of the ice layer, given by Equation ( 6), is discretized in finite elements.The FEM is used to calculate the average electric field Eavg, the residual resistance R(x), and the peak value of the leakage current Im determined by integrating the current density along a moving line located at equal distance of the arc roots.The FEM calculation are performed at each time step of the dynamic calculation algorithm, which used the numerical FEM results to verify the re-ignition condition given by Equation ( 16) followed by the arc propagation condition given by Equation (9).For each case, if the condition is not satisfied, the peak value of the applied source voltage Vm is increased by 1 kV, and the process is repeated until condition verification.Once the two conditions are verified, the velocity of each partial arc using Equation ( 10) is determined and used to increment the arc root position of these latter.The process is repeated until the sum of the length of the two partial arcs equals the arcing distance of the insulator.The simulation parameters used by the calculation algorithm are summarized in Table 1. Figure 6 shows the calculation algorithm of the dynamic numerical bi-arc models.This algorithm is based on the calculation algorithm developed for the dynamic numerical single arc presented in detail in [23] with the addition of the new conditions expressed by Equations ( 14)-( 16).The problem is solved in 2D using the electric current module and the stationary solver of Comsol Multiphysics, where only the conductive surface of the ice layer, given by Equation ( 6), is discretized in finite elements.The FEM is used to calculate the average electric field E avg , the residual resistance R(x), and the peak value of the leakage current I m determined by integrating the current density along a moving line located at equal distance of the arc roots.The FEM calculation are performed at each time step of the dynamic calculation algorithm, which used the numerical FEM results to verify the re-ignition condition given by Equation ( 16) followed by the arc propagation condition given by Equation ( 9).For each case, if the condition is not satisfied, the peak value of the applied source voltage V m is increased by 1 kV, and the process is repeated until condition verification.Once the two conditions are verified, the velocity of each partial arc using Equation ( 10) is determined and used to increment the arc root position of these latter.The process is repeated until the sum of the length of the two partial arcs equals the arcing distance of the insulator.The simulation parameters used by the calculation algorithm are summarized in Table 1.The model was validated using the experimental results extracted from literature [14,22].For that purpose, post-type insulators as the one in Figure 5a, with arcing distances of 1.39 and 2.02 m, were used.The corresponding FEM model is presented in Figure 5b.The experimental FOV of the post insulator was determined under wet-grown ice conditions during the melting period, the most likely conditions for flashover [4,5].The tests were conducted for two freezing water conductivities, 30 and 80 µs/cm and two ice thicknesses, 1.5 and 2 cm.The two air gaps were created artificially in order to control their length, fixed at 5 cm, which is the length used in the simulations.
Table 2 presents the comparison between the experimental results and those calculated with the proposed bi-arc numerical model, showing that the dynamic numerical model provides accurate results with an average discrepancy of 1.4%.Moreover, the implementation of the arc velocity criterion allows for computation of the evolution of the velocity of each partial arc during the flashover process presented in Figures 7 and 8 and the leakage current presented on Figure 9.The results were obtained for an arcing distance of 139 and 202 cm for an ice thickness of 1.5 cm and for freezing water conductivities of 30 and 80 µs/cm, respectively.As observed in Figures 7 and 8, only the velocity of one partial arc is represented.Indeed, Partial Arc 1 propagates downward and Partial Arc 2 propagates upward at the same velocity.Consequently, the flashover occurs when the two partial arcs meet each other at the middle of the ice layer or when they reach the same arc length equal to 64.5 cm for an arcing distance of 139 cm and equal to 96 cm for an arcing distance of 202 cm (Figure 8).The identical behavior of the two partial arcs can be attributed to the fact that the simulations were conducted for an identical initial arc length (x 01 = x 02 = 5 cm) for each partial arc.In this way, the problem becomes symmetric in terms of the arc velocity criterion given by Equation (10), as the two partial arcs are in series (with the same leakage current I m ) and each arc root has the same capacity C i calculated by Equation (12), between the arc root and the middle of the ice layer.Additionally, the velocity of each partial arc for the same arcing distance increases as freezing water conductivity increases.In the same way, the velocity of each partial arc, for the same freezing water conductivity, increases as arcing distance increases.From Figure 7, it can be observed that the time to flashover decreases as freezing water conductivity increases, for each arcing distance.The first observation of Figure 9 shows that the evolution and order of magnitude of the values obtained for the leakage current are consistent with the experimental observations made in laboratory [12][13][14][15][16][17][18].The results show that the leakage current increases with an increase in freezing water conductivity as well as with an increase in arcing distance.It can also be observed that the time to flashover decreases significantly as freezing water conductivity increases.The first observation of Figure 9 shows that the evolution and order of magnitude of the values obtained for the leakage current are consistent with the experimental observations made in laboratory [12][13][14][15][16][17][18].The results show that the leakage current increases with an increase in freezing water conductivity as well as with an increase in arcing distance.It can also be observed that the time to flashover decreases significantly as freezing water conductivity increases.The first observation of Figure 9 shows that the evolution and order of magnitude of the values obtained for the leakage current are consistent with the experimental observations made in laboratory [12][13][14][15][16][17][18].The results show that the leakage current increases with an increase in freezing water conductivity as well as with an increase in arcing distance.It can also be observed that the time to flashover decreases significantly as freezing water conductivity increases.

The Effect of the Initial Partial Arc Length
In order to verify the effect of the initial arc length, it was decided to perform another simulation, with the same arcing distances and parameters, but with the initial length of Partial Arc 1 equal to twice the length x2 of Partial Arc 2 (x01 = 2.x02).The arc velocity evolution of each partial arc as a function of time and the arc length for an arcing distance of 202 cm and freezing water conductivities of 30 and 80 µs/cm are presented in Figures 10 and 11, respectively.The results obtained clearly demonstrate that the initial arc length has a significant influence on its velocity.Indeed, it is interesting to note that the velocity of Arc 1 is significantly greater than that of Partial Arc 2, with a final value twice the velocity of Arc 2 for the arcing distance of 139 cm, the latter being similar to the velocity obtained with the same initial arc length (x01 = x02).For the arcing distance of 202 cm, the ratio between the velocities of Partial Arcs 1 and 2 decreases to a value of 1.3.With a higher velocity, Partial Arc 1 propagates on a longer distance at the ice surface than Partial Arc 2. In the case of Figure 10, Partial Arc 1 meets with Partial Arc 2 at a length x1 of 80 and 111 cm for the arcing distance of 139 and 202 cm, respectively (Figure 11).

The Effect of the Initial Partial Arc Length
In order to verify the effect of the initial arc length, it was decided to perform another simulation, with the same arcing distances and parameters, but with the initial length x 1 of Partial Arc 1 equal to twice the length x 2 of Partial Arc 2 (x 01 = 2x 02 ).The arc velocity evolution of each partial arc as a function of time and the arc length for an arcing distance of 202 cm and freezing water conductivities of 30 and 80 µs/cm are presented in Figures 10 and 11, respectively.The results obtained clearly demonstrate that the initial arc length has a significant influence on its velocity.Indeed, it is interesting to note that the velocity of Arc 1 is significantly greater than that of Partial Arc 2, with a final value twice the velocity of Arc 2 for the arcing distance of 139 cm, the latter being similar to the velocity obtained with the same initial arc length (x 01 = x 02 ).For the arcing distance of 202 cm, the ratio between the velocities of Partial Arcs 1 and 2 decreases to a value of 1.3.With a higher velocity, Partial Arc 1 propagates on a longer distance at the ice surface than Partial Arc 2. In the case of Figure 10, Partial Arc 1 meets with Partial Arc 2 at a length x 1 of 80 and 111 cm for the arcing distance of 139 and 202 cm, respectively (Figure 11).

The Effect of the Initial Partial Arc Length
In order to verify the effect of the initial arc length, it was decided to perform another simulation, with the same arcing distances and parameters, but with the initial length x1 of Partial Arc 1 equal to twice the length x2 of Partial Arc 2 (x01 = 2.x02).The arc velocity evolution of each partial arc as a function of time and the arc length for an arcing distance of 202 cm and freezing water conductivities of 30 and 80 µs/cm are presented in Figures 10 and 11, respectively.The results obtained clearly demonstrate that the initial arc length has a significant influence on its velocity.Indeed, it is interesting to note that the velocity of Arc 1 is significantly greater than that of Partial Arc 2, with a final value twice the velocity of Arc 2 for the arcing distance of 139 cm, the latter being similar to the velocity obtained with the same initial arc length (x01 = x02).For the arcing distance of 202 cm, the ratio between the velocities of Partial Arcs 1 and 2 decreases to a value of 1.3.With a higher velocity, Partial Arc 1 propagates on a longer distance at the ice surface than Partial Arc 2. In the case of Figure 10, Partial Arc 1 meets with Partial Arc 2 at a length x1 of 80 and 111 cm for the arcing distance of 139 and 202 cm, respectively (Figure 11).Table 3 presents a comparison between the FOV results obtained with the bi-arc dynamic numerical model for the two partial arc initial length conditions.As can be observed, the effect of initial arc length is not significant on the FOV results.However, the absolute discrepancy seems to increase with an increase of arcing distance and freezing water conductivity as well.

Implementation of the Mayr Approach in the Bi-Arc Dynamic Numerical Model
As shown in the last section, the proposed numerical dynamic model improves the accuracy of the initial numerical static bi-arc model and computes the evolution of the leakage current and the arc velocity during the flashover process.However, as explained in Section 3.2, the proposed calculation algorithm is applicable only in the case where the two partial arcs are in contact with the HV and ground electrodes.In this particular situation, the voltage condition applied to each arc root can be easily determined using Equations ( 14) and (15).In the case of one partial arc having its two arc roots in contact with the ice layer, as illustrated by Arc 2 in Figure 4a, these equations are no longer valid.The two arc roots must be considered as two floating potentials of unknown values, which are required to deeply modify the calculation algorithm of Figure 6.In order to solve this problem, it was decided to implement a new modeling of the arc in the numerical dynamic bi-arc model using the arc resistance formulation proposed by Mayr, as presented in detail in the next section.Table 3 presents a comparison between the FOV results obtained with the bi-arc dynamic numerical model for the two partial arc initial length conditions.As can be observed, the effect of initial arc length is not significant on the FOV results.However, the absolute discrepancy seems to increase with an increase of arcing distance and freezing water conductivity as well.

Implementation of the Mayr Approach in the Bi-Arc Dynamic Numerical Model
As shown in the last section, the proposed numerical dynamic model improves the accuracy of the initial numerical static bi-arc model and computes the evolution of the leakage current and the arc velocity during the flashover process.However, as explained in Section 3.2, the proposed calculation algorithm is applicable only in the case where the two partial arcs are in contact with the HV and ground electrodes.In this particular situation, the voltage condition applied to each arc root can be easily determined using Equations ( 14) and (15).In the case of one partial arc having its two arc roots in contact with the ice layer, as illustrated by Arc 2 in Figure 4a, these equations are no longer valid.The two arc roots must be considered as two floating potentials of unknown values, which are required to deeply modify the calculation algorithm of Figure 6.In order to solve this problem, it was decided to implement a new modeling of the arc in the numerical dynamic bi-arc model using the arc resistance formulation proposed by Mayr, as presented in detail in the next section.Mayr has proposed an analytical formulation in order to define the arc resistance R arc , which can be defined only by the arc voltage and current.He has also described the arc resistance as time-dependent by introducing the arc time constant τ.According to Mayr, the arc resistance R arc can be expressed as follows [26]: where τ = 100 µs is the arc time constant [26] and P 0 the cooling power dissipated in the arc, expressed numerically as where I m is the leakage current (A), x is the arc length, and E arc is the electric field along the arc given by Equation ( 2).The development of Equation ( 17) for its implementation in the calculation algorithm leads to the following expression: In order to be consistent with the implementation of the arc resistance formulation given by Equation (19) in the algorithm, the authors decided to introduce a new expression of the electric field E arc along the partial arc, which is generally expressed by the Nottingham expression given by Equation (2).The new expression of E arc is as follows: where I m is le leakage current, R arc is the arc resistance given by Equation ( 19), and x is the partial arc length.

Validation of the Implementation of the Mayr Arc Resistance Formulation
The implementation of Equations ( 19) and (20) in the dynamic bi-arc algorithm of Figure 6 was validated using the same simulation parameters and ice-covered insulator geometries as used in Section 3.3.The results obtained were then compared with those of the Nottingham approach in terms of FOV, arc velocity, and leakage current evolutions.Table 4 presents the FOV predictive results obtained with the two approaches for an ice thickness of 1.5 cm, which are compared to the experimental results.The results showed that the Mayr approach is less accurate than Nottingham with an average discrepancy of 2.43% and 0.86%, respectively.However, the difference between the two approaches is not significant.The absolute accuracy seems to decrease significantly as arcing distance and freezing water conductivity are increased, but less so for the latter.
Figures 12 and 13 compare the evolution of partial arc velocity as a function of time and arc length, respectively, for an arcing distance of 202 cm and initial condition x 01 = x 02 .Only the velocity of one partial arc is presented, as they are the same due to the symmetry of the problem, as explained previously.Figure 14 shows the evolution of leakage current I m as a function of arc length for an arcing distance of 202 cm and initial condition x 01 = x 02 .
From Figure 12, it can be observed that the evolution of the velocities obtained with Mayr is different from those obtained with Nottingham, with a longer time to flashover for Mayr.In Figure 13, however, the evolutions of the arc velocity as a function of the arc length are closed, as the dynamic aspect is not represented in this figure, with a final value slightly greater for Nottingham.The same observation goes for the evolution of leakage current in Figure 14 but with a difference in the final value of I m , which is higher for Mayr than it is for Nottingham.Finally, the authors decided to compare the effect of the initial arc length between the Mayr and Nottingham approaches.FOV results are listed in Table 5 for an ice thickness of 1.5 cm, an arcing distance of 202 cm, and a freezing water conductivity of 80 µs/cm.With an average discrepancy of 2.37%, compared to the 1.53% obtained with Nottingham, it can be concluded that the Mayr approach, in the same manner that Nottingham, takes into account the effect of the initial partial arc length.
Figures 15 and 16 present the evolution of the velocity of each partial arc as a function of time and arc length, respectively.The simulations were performed for the same simulation conditions used for an arcing distance of 202 cm and for a freezing water conductivity of 80 µs/cm.As observed previously, the dynamic behavior of the partial arc differs from Mayr and Nottingham with a greater time to flashover with the Mayr approach.However, the arc velocity values obtained with Mayr are in good agreement with those obtained with Nottingham, which is confirmed by the results of Figure 16.Finally, the authors decided to compare the effect of the initial arc length between the Mayr and Nottingham approaches.FOV results are listed in Table 5 for an ice thickness of 1.5 cm, an arcing distance of 202 cm, and a freezing water conductivity of 80 µs/cm.With an average discrepancy of 2.37%, compared to the 1.53% obtained with Nottingham, it can be concluded that the Mayr approach, in the same manner that Nottingham, takes into account the effect of the initial partial arc length.
Figures 15 and 16 present the evolution of the velocity of each partial arc as a function of time and arc length, respectively.The simulations were performed for the same simulation conditions used for an arcing distance of 202 cm and for a freezing water conductivity of 80 µs/cm.As observed previously, the dynamic behavior of the partial arc differs from Mayr and Nottingham with a greater time to flashover with the Mayr approach.However, the arc velocity values obtained with Mayr are in good agreement with those obtained with Nottingham, which is confirmed by the results of Figure 16.

Validation of the Bi-Arc Dynamic Numerical Model Results
The obtained results have demonstrated the predicting ability of the proposed bi-arc dynamic numerical model with an average discrepancy of 1.4% for the FOV values (Table 2).In that case, the initial length of the two air gaps were the same and equal to 5 cm in order to respect the experimental conditions used in the literature to determine the FOV of ice-covered post-type insulators.This result can also be compared to that obtained with the authors' previous bi-arc static numerical model [22] with which an average discrepancy of 4.2% was obtained for the same experimental conditions.This clearly demonstrates that the implementation of the Galimberti arc velocity criterion improves the accuracy of the FOV predictive results.Moreover, this implementation also allows one to compute the velocity of each partial arc during the flashover

Validation of the Bi-Arc Dynamic Numerical Model Results
The obtained results have demonstrated the predicting ability of the proposed bi-arc dynamic numerical model with an average discrepancy of 1.4% for the FOV values (Table 2).In that case, the initial length of the two air gaps were the same and equal to 5 cm in order to respect the experimental conditions used in the literature to determine the FOV of ice-covered post-type insulators.This result can also be compared to that obtained with the authors' previous bi-arc static numerical model [22] with which an average discrepancy of 4.2% was obtained for the same experimental conditions.This clearly demonstrates that the implementation of the Galimberti arc velocity criterion improves the accuracy of the FOV predictive results.Moreover, this implementation also allows one to compute the velocity of each partial arc during the flashover

Validation of the Bi-Arc Dynamic Numerical Model Results
The obtained results have demonstrated the predicting ability of the proposed bi-arc dynamic numerical model with an average discrepancy of 1.4% for the FOV values (Table 2).In that case, the initial length of the two air gaps were the same and equal to 5 cm in order to respect the experimental conditions used in the literature to determine the FOV of ice-covered post-type insulators.This result can also be compared to that obtained with the authors' previous bi-arc static numerical model [22] with which an average discrepancy of 4.2% was obtained for the same experimental conditions.This clearly demonstrates that the implementation of the Galimberti arc velocity criterion improves the accuracy of the FOV predictive results.Moreover, this implementation also allows one to compute the velocity of each partial arc during the flashover process, as shown in Figures 7 and 8. Due to the symmetrical configuration of the problem when the two partial arcs present the same initial length, the two partial arcs propagate with the same velocity.This velocity is influenced by the freezing water conductivity by around 30%.The most significant difference between the two approaches lies in the simulation of the dynamic behavior of the partial arcs in terms of time to flashover, which is longer in the case of Mayr, as can be observed in Figures 12 and 15.This difference can be attributed to the expression of the arc resistance which remains an approximation of the arc behavior, and which was initially developed for polluted insulators [26].Moreover, the increase in time to flashover obtained with Mayr can be explained by the implementation of the expression of E arc by Equation (20), which is dependent on the arc resistance calculation.However, despite the difference obtained, the results obtained with the Mayr approach can be considered accurate, taking into account the effect of freezing water conductivity in the same manner as the Nottingham approach.

The Proposition to Extend the Bi-Arc Dynamic Numerical Model to Multi-Arc Flashover Modeling
In this section, the possibility to extend the bi-arc dynamic numerical model to a multi-arc flashover process is discussed.As illustrated in Figure 3, the flashover of EHV post-type ice-covered insulators with an arcing distance of up to 4.17 m can be the result of the propagation of three to several partial arcs.In general, when the arcing distance becomes greater than 2.02 m, at least one of the arcs has its two arc roots in contact with the ice layer, as illustrated by Arc 2 in Figure 17a.This arc configuration is problematic, as its two arc roots become two equipotential surfaces with floating potentials, adding two unknown values in the electrical circuit.However, this specific constraint can be solved by using the Mayr approach, as given by Equation ( 19), to model Partial Arc 2 as an electrical resistance.
sensibly less accurate than those obtained by Nottingham with an average discrepancy of 2.4% compared to the 1.2% obtained with the latter.In Figures 12-16, it can be observed that Mayr seems to underestimate the final value of the arc velocity by around 7% but overestimates the final value of the leakage current by around 30%.The most significant difference between the two approaches lies in the simulation of the dynamic behavior of the partial arcs in terms of time to flashover, which is longer in the case of Mayr, as can be observed in Figures 12 and 15.This difference can be attributed to the expression of the arc resistance which remains an approximation of the arc behavior, and which was initially developed for polluted insulators [26].Moreover, the increase in time to flashover obtained with Mayr can be explained by the implementation of the expression of Earc by Equation (20), which is dependent on the arc resistance calculation.However, despite the difference obtained, the results obtained with the Mayr approach can be considered accurate, taking into account the effect of freezing water conductivity in the same manner as the Nottingham approach.

The Proposition to Extend the Bi-Arc Dynamic Numerical Model to Multi-Arc Flashover Modeling
In this section, the possibility to extend the bi-arc dynamic numerical model to a multi-arc flashover process is discussed.As illustrated in Figure 3, the flashover of EHV post-type ice-covered insulators with an arcing distance of up to 4.17 m can be the result of the propagation of three to several partial arcs.In general, when the arcing distance becomes greater than 2.02 m, at least one of the arcs has its two arc roots in contact with the ice layer, as illustrated by Arc2 in Figure 17a.This arc configuration is problematic, as its two arc roots become two equipotential surfaces with floating potentials, adding two unknown values in the electrical circuit.However, this specific constraint can be solved by using the Mayr approach, as given by Equation ( 19), to model Partial Arc 2 as an electrical resistance.By using the Mayr approach for Partial Arc 2, the multi-arc model can be divided into several bi-arc modules as illustrated in Figure 17b.The bi-arc calculation algorithm with the implementation of an arc resistance calculation can then be applied to each module.Each bi-arcs module is in series and interconnected through the partial arc with its two arc roots in contact with the ice layer (Arc 2 in Figure 17).Hence, the modules are interdependent; consequently, they have an influence on each other.Figure 18 proposes a calculation algorithm that can be used to develop a new multi-arc dynamic numerical model involving three partial arcs at the ice surface, as shown in Figure 17a.In this example, each bi-arc module calculates the different parameters associated to their respective ice layer section (arc propagation, arc velocity, arc and residual resistances, and Eavg) and exchanges these data to compute the leakage current Im as a function of applied voltage Vm.Once this is done, the algorithm verifies if one of the sections of the ice layer is short-circuited by the two partial arcs (partial flashover obtained).In the negative, the process is repeated until a partial flashover occurs.When this happens, the multi-arc algorithm applies the bi-arc module to the remaining ice section By using the Mayr approach for Partial Arc 2, the multi-arc model can be divided into several bi-arc modules as illustrated in Figure 17b.The bi-arc calculation algorithm with the implementation of an arc resistance calculation can then be applied to each module.Each bi-arcs module is in series and interconnected through the partial arc with its two arc roots in contact with the ice layer (Arc 2 in Figure 17).Hence, the modules are interdependent; consequently, they have an influence on each other.Figure 18 proposes a calculation algorithm that can be used to develop a new multi-arc dynamic numerical model involving three partial arcs at the ice surface, as shown in Figure 17a.In this example, each bi-arc module calculates the different parameters associated to their respective ice layer section (arc propagation, arc velocity, arc and residual resistances, and E avg ) and exchanges these data to compute the leakage current I m as a function of applied voltage V m .Once this is done, the algorithm verifies if one of the sections of the ice layer is short-circuited by the two partial arcs (partial flashover obtained).In the negative, the process is repeated until a partial flashover occurs.When this happens, the multi-arc algorithm applies the bi-arc module to the remaining ice section with the two remaining partial arcs.In a general way, when N partial arcs are present at the ice surface with N−2 partial arcs having their two arc roots in contact with the ice layer, N−1 bi-arc modules are required to simulate this problem using the same principle of calculation described previously.In this way, the proposed calculation algorithm of Figure 18 has no limit in terms of the number of partial arcs.
with the two remaining partial arcs.In a general way, when N partial arcs are present at the ice surface with N−2 partial arcs having their two arc roots in contact with the ice layer, N−1 bi-arc modules are required to simulate this problem using the same principle of calculation described previously.In this way, the proposed calculation algorithm of Figure 18 has no limit in terms of the number of partial arcs.It is important to mention that, as illustrated in Figure 17b, each bi-arc module takes into account the entire length of the partial arc (Arc 2 in the example of Figure 17) in the validation of the arc propagation and arc velocity criteria given by Equations ( 9) and (10).The length of Partial Arc 2 is important because, as demonstrated by the results of Figures 10 and 11, it has a significant influence on its velocity, in particular.This approach differs from that used in the multi-arc dynamic mathematical model and illustrated in Figure 4 where Partial Arc 2 is divided into two partial arcs of equal length.However, this mathematical model makes use of an experimental arc velocity criterion given by Equation (8), which only allows one to determine the average velocity of the partial arc and not the velocity of each partial arc, as reported in [18].Moreover, this mathematical model uses the modified Wilkins formulation given by Equation (7), which can only be applied to rectangular geometries, hence limiting considerably its applicability to more complex problems such as polluted insulators.

Conclusions
This paper proposes a new bi-arc dynamic numerical model based on the use of commercial finite element software coupled with a generic calculation algorithm.This new model can simulate the final stage of the flashover process of ice-covered insulators resulting from the propagation of two partial arcs at the ice surface.The accuracy of the results obtained, with an average discrepancy of 1.4% for the prediction of FOV, conclusively shows that the implementation of the Hampton and Gallimberti criteria in a modular calculation algorithm can obtain a reliable predictive model.It is important to mention that, as illustrated in Figure 17b, each bi-arc module takes into account the entire length of the partial arc (Arc 2 in the example of Figure 17) in the validation of the arc propagation and arc velocity criteria given by Equations ( 9) and (10).The length of Partial Arc 2 is important because, as demonstrated by the results of Figures 10 and 11, it has a significant influence on its velocity, in particular.This approach differs from that used in the multi-arc dynamic mathematical model and illustrated in Figure 4 where Partial Arc 2 is divided into two partial arcs of equal length.However, this mathematical model makes use of an experimental arc velocity criterion given by Equation (8), which only allows one to determine the average velocity of the partial arc and not the velocity of each partial arc, as reported in [18].Moreover, this mathematical model uses the modified Wilkins formulation given by Equation (7), which can only be applied to rectangular geometries, hence limiting considerably its applicability to more complex problems such as polluted insulators.

Conclusions
This paper proposes a new bi-arc dynamic numerical model based on the use of commercial finite element software coupled with a generic calculation algorithm.This new model can simulate the final stage of the flashover process of ice-covered insulators resulting from the propagation of two partial arcs at the ice surface.The accuracy of the results obtained, with an average discrepancy of 1.4% for the prediction of FOV, conclusively shows that the implementation of the Hampton and Gallimberti criteria in a modular calculation algorithm can obtain a reliable predictive model.Moreover, thanks to the use of the finite element method, it is now possible to avoid the problem of analytical formulation for the residual resistance calculation, consequently allowing the application of this new bi-arc model

Figure 1 .
Figure 1.Single arc modeling of flashover proposed by Obenaus.

Figure 1 .
Figure 1.Single arc modeling of flashover proposed by Obenaus.

Figure 2 .
Figure 2. Example of ice accumulation shape obtained for (a) a porcelain post insulator with an arcing distance L of 103 cm and (b) the corresponding rectangular model geometry used in the flashover model.

Figure 2 .
Figure 2. Example of ice accumulation shape obtained for (a) a porcelain post insulator with an arcing distance L of 103 cm and (b) the corresponding rectangular model geometry used in the flashover model.

Figure 3 .
Figure 3. Example of ice accumulation shape obtained for an extra-high voltage (EHV) porcelain post insulator under wet grown ice conditions: (a) with the presence of two partial arcs for an arcing distance of 1.36 m; (b) with the presence of three partial arcs for an arcing distance of 2.7 m [14,18].

Figure 4 .
Figure 4. Illustration of the principle use for modelling multi-arc flashover process of ice-covered EHV insulator: (a) initial model with three partial arcs; (b) equivalent model composed of four single arcs.

Figure 3 .
Figure 3. Example of ice accumulation shape obtained for an extra-high voltage (EHV) porcelain post insulator under wet grown ice conditions: (a) with the presence of two partial arcs for an arcing distance of 1.36 m; (b) with the presence of three partial arcs for an arcing distance of 2.7 m [14,18].

Figure 3 .
Figure 3. Example of ice accumulation shape obtained for an extra-high voltage (EHV) porcelain post insulator under wet grown ice conditions: (a) with the presence of two partial arcs for an arcing distance of 1.36 m; (b) with the presence of three partial arcs for an arcing distance of 2.7 m [14,18].

Figure 4 .
Figure 4. Illustration of the principle use for modelling multi-arc flashover process of ice-covered EHV insulator: (a) initial model with three partial arcs; (b) equivalent model composed of four single arcs.

Figure 4 .
Figure 4. Illustration of the principle use for modelling multi-arc flashover process of ice-covered EHV insulator: (a) initial model with three partial arcs; (b) equivalent model composed of four single arcs.

Figure 5 .
Figure 5. Modelling of the bi-arc flashover process of the ice-covered EHV insulator.

Figure 5 .
Figure 5. Modelling of the bi-arc flashover process of the ice-covered EHV insulator.

Figure 6 .
Figure 6.Calculation algorithm of the bi-arc dynamic numerical model.Figure 6. Calculation algorithm of the bi-arc dynamic numerical model.

Figure 6 .
Figure 6.Calculation algorithm of the bi-arc dynamic numerical model.Figure 6. Calculation algorithm of the bi-arc dynamic numerical model.

3. 3 .
The Validation of the Bi-Arc Dynamic Numerical Model3.3.1.The Effect of Arcing Distance and Freezing Water Conductivity

Figure 7 .
Figure 7. Evolution of partial arc velocity as a function of time for equal initial arc length (x01 = x02) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 8 .
Figure 8. Evolution of partial arc velocity as a function of arc length for an equal initial arc length (x01 = x02) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 7 .
Figure 7. Evolution of partial arc velocity as a function of time for equal initial arc length (x 01 = x 02 ) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 7 .
Figure 7. Evolution of partial arc velocity as a function of time for equal initial arc length (x01 = x02) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 8 .
Figure 8. Evolution of partial arc velocity as a function of arc length for an equal initial arc length (x01 = x02) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 8 .
Figure 8. Evolution of partial arc velocity as a function of arc length for an equal initial arc length (x 01 = x 02 ) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 9 .
Figure 9. Evolution of the leakage current Im as a function of time for an equal initial arc length (x01 = x02) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 10 .
Figure 10.Evolution of partial arc velocities as a function of arc length for two different initial arc lengths (x01 = 2x02) obtained for an arcing distance of 139 cm with a freezing water conductivity of 80 µs/cm.

Figure 9 .
Figure 9. Evolution of the leakage current I m as a function of time for an equal initial arc length (x 01 = x 02 ) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Energies 2016, 9 , 22 Figure 9 .
Figure 9. Evolution of the leakage current Im as a function of time for an equal initial arc length (x01 = x02) obtained for arcing distances of 139 and 202 cm with freezing water conductivities of 30 and 80 µs/cm.

Figure 10 .
Figure 10.Evolution of partial arc velocities as a function of arc length for two different initial arc lengths (x01 = 2x02) obtained for an arcing distance of 139 cm with a freezing water conductivity of 80 µs/cm.

Figure 10 .
Figure 10.Evolution of partial arc velocities as a function of arc length for two different initial arc lengths (x 01 = 2x 02 ) obtained for an arcing distance of 139 cm with a freezing water conductivity of 80 µs/cm.

Figure 11 .
Figure 11.Evolution of the partial arc velocities as a function of arc length for two different initial arc lengths (x01 = 2x02) obtained for an arcing distance of 202 cm with a freezing water conductivity of 80 µs/cm.

Figure 11 .
Figure 11.Evolution of the partial arc velocities as a function of arc length for two different initial arc lengths (x 01 = 2x 02 ) obtained for an arcing distance of 202 cm with a freezing water conductivity of 80 µs/cm.

1 .
Arc Resistance Formulation by Mayr

Figure 12 .
Figure 12.Evolution of partial arc velocities as a function of time between the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x01 = x02).

Figure 13 .
Figure 13.Evolution of partial arc velocities as a function of the arc length between the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x01 = x02).

Figure 12 .
Figure 12.Evolution of partial arc velocities as a function of time between the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x 01 = x 02 ).

Figure 12 .
Figure 12.Evolution of partial arc velocities as a function of time between the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x01 = x02).

Figure 13 .
Figure 13.Evolution of partial arc velocities as a function of the arc length between the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x01 = x02).

Figure 13 .
Figure 13.Evolution of partial arc velocities as a function of the arc length between the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x 01 = x 02 ).

Figure 14 .
Figure 14.Evolution of leakage current Im as a function of arc length obtained with the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x01 = x02).

Figure 14 .
Figure 14.Evolution of leakage current I m as a function of arc length obtained with the Nottingham and Mayr approaches for an arcing distance of 202 cm and the same initial arc length (x 01 = x 02 ).

Figure 15 .
Figure 15.Evolution of the each arc velocity as a function of time with the Nottingham and Mayr approaches for an arcing distance of 202 cm, a freezing water conductivity of 80 µs/cm, and an initial arc length x01 = 2x02.

Figure 16 .
Figure 16.Evolution of the each arc velocity as a function of arc length with the Nottingham and Mayr approaches for an arcing distance of 202 cm, a freezing water conductivity of 80 µs/cm, and an initial arc length x01 = 2x02.

Figure 15 . 22 Figure 15 .
Figure 15.Evolution of the each arc velocity as a function of time with the Nottingham and Mayr approaches for an arcing distance of 202 cm, a freezing water conductivity of 80 µs/cm, and an initial arc length x 01 = 2x 02 .

Figure 16 .
Figure 16.Evolution of the each arc velocity as a function of arc length with the Nottingham and Mayr approaches for an arcing distance of 202 cm, a freezing water conductivity of 80 µs/cm, and an initial arc length x01 = 2x02.

Figure 16 .
Figure 16.Evolution of the each arc velocity as a function of arc length with the Nottingham and Mayr approaches for an arcing distance of 202 cm, a freezing water conductivity of 80 µs/cm, and an initial arc length x 01 = 2x 02 .

Figure 17 .
Figure 17.Examples of decomposition of a multi-arc problem with several bi-arc modules: (a) multi-arc models of arcing distance L with three partial arcs and (b) decomposition of the model in two bi-arc modules.

Figure 17 .
Figure 17.Examples of decomposition of a multi-arc problem with several bi-arc modules: (a) multi-arc models of arcing distance L with three partial arcs and (b) decomposition of the model in two bi-arc modules.

Figure 18 .
Figure 18.Principle of the calculation algorithm developed for the numerical dynamic multi-arc models in the case of three partial arcs present at the ice surface.

Figure 18 .
Figure 18.Principle of the calculation algorithm developed for the numerical dynamic multi-arc models in the case of three partial arcs present at the ice surface.

Table 2 .
Flashover voltage (FOV) obtained with the new numerical dynamic bi-arc model for a freezing water conductivity of 80 µs/cm.

Table 3 .
Effect of initial arc length on the FOV predictive results for freezing water conductivity of 80 µs/cm.

Table 3 .
Effect of initial arc length on the FOV predictive results for freezing water conductivity of 80 µs/cm.

Table 4 .
FOV predictive results obtained with the Nottingham and Mayr approaches for an ice thickness of 1.5 cm and the same initial arc length (x 01 = x 02 ).

Table 4 .
FOV predictive results obtained with the Nottingham and Mayr approaches for an ice thickness of 1.5 cm and the same initial arc length (x01 = x02).

Table 4 .
FOV predictive results obtained with the Nottingham and Mayr approaches for an ice thickness of 1.5 cm and the same initial arc length (x01 = x02).

Table 5 .
FOV predictive results obtained with the Nottingham and Mayr approaches for an ice thickness of 1.5 cm and for different initial arc lengths (x01 = 2x02).

Table 5 .
FOV predictive results obtained with the Nottingham and Mayr approaches for an ice thickness of 1.5 cm and for different initial arc lengths (x 01 = 2x 02 ).