High Voltage Special Issue: Interface Charging Phenomena for Dielectric Materials Finite-element analysis for surface discharge on polyimide insulation in air at atmospheric pressure under pulsed electrical stress

A surface discharge non-equilibrium plasma model of air–polyimide under pulsed electrical stress is established, by considering the reaction of charged particles on the dielectric surface and the secondary electron emission caused by the condition that high-energy particles bombard the material surface. The model defines the chemical reaction of air discharge by using simplified set of reactions, which greatly reduces the complexity of the model. To avoid the negative value of particle density in the process of solution, the logarithmic finite-element method is used to solve the model established, so as to implement the dynamic simulation of the surface discharge process. Also, the temporal and spatial evolution of the microparameters such as charge and electric field distribution during discharge are obtained, and the reliability of the model is verified by experiments in terms of discharge development pattern and surface charge accumulation. By comparing the development process of surface discharge under single pulse and repetitive pulses, it can be seen that surface discharge develops from needle electrode to ground electrode under both repetitive pulses and single pulse stress, but the relationship between the discharge propagation time under repetitive pulses and pulse repetition rate is a ‘u’ curve, and the inflection point moves to higher repetition rate region with the increase of voltage.


Introduction
As a new type of electric transmission and transformation equipment with flexible transformation mode and excellent control performance, high-frequency power transformer (HFPT) is of promising application prospects in electric traction, distributed energy grid-connected and other fields [1]. HFPT mainly uses airpolyimide (PI) insulation system [2]. Due to its small size, compact structure, long-term high frequency electrical stress with short rising time and large amplitude [3], the operating environment is worse than that of traditional transformers. This is easier to induce partial discharge or surface discharge at weak part of insulation, which results in ageing or premature failure of PI. Besides, with the increase of transformer capacity and voltage class, problems about insulation become more prominent [4,5].
The stable operation of high-frequency power equipment largely depends on its insulation level, and the gas-solid interface is the weak part of insulation. With the development of EHV/UHV power transmission, surface discharge has become one of the important reasons that lead to the failure of gas-solid insulation system of high-frequency equipment and endanger the security of power system [6,7]. Due to the limitation of existing experimental methods, the distribution of charged particles and electric field distortion during surface discharge are difficult to measure, and the micro-mechanism of discharge process cannot be fully clarified [8]. With the leap of computer computing and storage capacity and development of numerical simulation methods, it is possible to simulate micro-parameters, such as electric field and charge distribution, in discharge process. At present, the simulation of discharge mostly focuses on air gap discharge [9,10] between needle-plate electrodes and dielectric barrier discharge (DBD) between plate electrodes [11,12]. Compared with surface discharge, the effect of solid insulation on the discharge process need not be considered in the air gap discharge process. DBD aims to produce a large area of low-temperature plasma with high energy density under certain conditions, and it will not form a penetrating discharge channel between the two electrodes. Currently, the simulation of the surface discharge process is still in the exploratory stage. Ho-Young Lee et al. [13] simulated the surface discharge process of the liquid-solid insulation system with needle-bar electrode structure by using the electro-hydrodynamics method and calculated the discharge propagation velocity. Based on the hydrodynamic theory, Shailendra Singh et al. [14] considered the charge transport process inside a dielectric barrier and analysed the influence of solid dielectric barrier on development of air discharge between needle and plate electrodes. Some scholars combined hydrodynamic theory with plasma reaction in discharge. For example, Wenxia Sima et al. [15] considered 12 types of particles and 25 types of chemical reactions to establish the surface discharge plasma model in N 2 /O 2 environment, so as to obtain distribution law of surface charge and other micro-parameters. Some progress has been made in these studies, but the electrode structure is mostly confined to the surface DBD (linear or comb electrodes distributed on both sides of dielectric plate) and the plate electrodes, while the surface discharge of high frequency power equipment mostly occurs in the extremely uneven electric field, and the extensive use of switching devices makes the insulation withstand repeated pulse stresses for long periods of time. It is necessary to carry out the simulation study of discharge surface of needle-bar electrode structure under pulse stress. In addition, the complex air components lead to hundreds of particle reactions in the discharge process, and how to rationally simplify the reactions still needs to be further explored.
In order to investigate the evolution process of surface discharge, a simplified but effective set of reactions is used to describe all the particle reactions in air discharges. On this basis, a numerical model for discharge along air-PI surface is proposed with fluid dynamic theory. Then, the surface discharge under pulsed electric stress is investigated, and the distribution of microparameters such as surface charge and particle density are obtained. In addition, the differences and similarities of surface discharge process under single pulse and repetitive pulses are analysed. All of these are helpful to reveal the micro-physical process of surface discharge.

Electrode structure
Nowadays, needle-plate [16] and needle-bar [17] are two main electrode structures used to simulate surface discharge. The strong normal component of electric field with the needle-plate electrode structure can easily cause partial discharge in the dielectric, leading to the longitudinal breakdown of solid insulation. The needle-bar model can change the tangential and normal electric field distribution by adjusting the inclination angle β of needle electrode, so as to flexibly control the propagation velocity of surface discharge, and facilitate the observation of discharge process. Therefore, the needle-bar electrode structure shown in Fig. 1a is chosen in this paper to simulate surface discharge behaviour. It is tested that when the angle between the needle electrode and PI film surface is 45°, there are sufficient tangential electric field components along the insulation surface and the normal electric field will not lead to thermal breakdown of insulation. The volume and weight of HFPT will be significantly reduced due to the increase of frequency. At present, most experiments to simulate the surface discharge of HFPT are based on 1-2 cm of gap [13,18]. The influence of surface discharge distance on the discharge inception voltage and propagation velocity are considered comprehensively, and the discharge gap d sg is set to 1 cm between needle and the bar electrode. Experimental results show that the discharge trajectory along PI under needle-bar electrode structure is approximately linear, in other words, the plasma channel mainly develops along the vertical line between the needle tip and the ground electrode [19]. To reduce the complexity of the model, the two-dimensional model shown in Fig. 1b is used to approximate the needle-bar electrode structure. The diameter of needle electrode is 2 mm, the radius of curvature of tip is 6 μm and the specifications of bar electrode are: length × width × height = 50 mm × 5 mm × 6 mm. The thickness of PI film d d = 125 μm, relative dielectric constant ε 0 = 1, ε rd = 3.7.

Governing equations
Surface discharge process not only needs to describe the discharge process in air accurately, but also needs to consider the influence of solid insulation. The traditional fluid model involves continuity equation of electron, positive and negative ions with Poisson equation, describing the collision, diffusion and drift of the charged particles under electric field and the self-consistent field produced by the charged particle. However, the secondary electron emission process caused by high-energy particles bombarding the material surface is not considered in the traditional fluid model. Therefore, it is necessary to add the secondary electron emission on the insulation surface as a electron source term on the basis of the traditional fluid equation.
In general, electron transport is described by the Boltzmann equation, which is a six-dimensional, non-local continuity equation in phase space. The Boltzmann equation is an extremely complicated integrodifferential equation and solving it in an efficient manner is not currently possible. However, the Boltzmann equation can be approximated by fluid equations by multiplying by a weighting function and then integrating over velocity space. This reduces the governing equations to a three-dimensional, timedependent problem.
The electron transport equation is ∂ ∂t n e + ∇ ⋅ Γ e = R e + S e (1) where n e is the electron density; Γ e is the electron flux vector; R e is either a source or a sink of electrons which is determined by plasma reaction and dielectric surface reaction; μ e is the electron mobility; D e is the electron diffusivity, E is electric field; S e is the photoelectron generation rate, and the traditional integral model can be approximately replaced by the three-exponential Helmholtz equations. The specific expressions and parameter values are derived from literature [20]. The rate of change of the electron energy density is described by where n ɛ is the electron energy density; S en is the energy loss or gain due to inelastic collisions; μ ɛ and D ɛ are, respectively, the electron energy mobility and the electron energy diffusivity. According to Einstein formula, the electron mobility, the electron diffusivity, the electron energy mobility and electron energy diffusivity satisfy the relationships: D e = μ e T e , μ ɛ = 5/3μ e and D ɛ = μ ɛ T e . Where, μ e = 3.74 × 10 24 (10 21 E) −0.22 (V*m*s) −1 [21]; T e is the electron temperature and the relationship between the mean electron energy ε ave and the electron temperature T e is T e = 2/3ε ave .
The source terms R e and S en in (1) and (3) are determined by plasma chemical reactions. Suppose that there are M(M = M 1 + M 2 ) reactions which contribute to the growth or decay of electron density and P(P = P 1 + P 2 ) inelastic electron-neutral collisions. When the reaction rate or Townsend coefficient is given, R e and S en can be expressed as where x i and x j are, respectively, the mole fractions of the target species for reaction i and reaction j; k i is the rate coefficient for reaction i; α j is the Townsend coefficient for reaction j; N n is the total neutral number density, and Δɛ i and Δɛ j are, respectively, the energy losses for reaction i and reaction j. Non-electron species can be described by the multi-component diffusion equation where k = 1, …, Q represents the type of species; ρ denotes the density of the mixture; w k is the mass fraction of the species k; u is the mass averaged fluid velocity vector; R k is the rate expression for species k; j k is the diffusive flux vector; S p is the photoionisation item, which only exists in the diffusion equation of According to the Poisson equation, the electric field could be calculated where ρ v is the space charge density; q represents the elementary charge; Z k is the quantity of the charge which k takes; n k is the density of k; V is the electric potential.

Chemical reactions
As a mixed gas, the discharge of air involves complicated plasma chemical reactions, including dozens kinds of particles and hundreds of chemical reactions [22,23], and it is almost impossible to establish a complete plasma reaction system. In fluid dynamics theory, various particles are described by positive, negative and neutral particles [9]. Considering the main goal of this model is to study the dynamic evolvement process of the discharge, and it is not concerned with the specific situation of a certain reaction, a simplified set of reactions is used to describe correctly the creation and destruction of charged species in air. The plasma reactions considered in this model are shown in Table 1. Nitrogen and oxygen are not treated separately. Instead a general species A is used for air. Species A can be ionised forming positive ions p and can also attach electrons forming negative ions n, where NA is the Avogadro constant. The unit of three-body reaction (reaction between three particles) rate coefficient is m 6 /s and that of two-body reaction rate coefficient is m 3 /s. The rates of ionisation and attachment reactions are given by Townsend coefficient and three-body attachment reaction, electron-positive ion and positive-negative ion recombination are determined by the rate coefficient [12]. Simplify the air components as N 2 :O 2 = 4:1. With the cross-section data of nitrogen and oxygen [24], Townsend coefficient is obtained by using Boltzmann two-term approximation method [25]. Considering only the primary ionisation of nitrogen and oxygen molecules, the rate coefficient k is where the constant γ = (2q e /m e ) 1/2 , and q e is the electron charge, m e is the electron mass and ε is energy; σ(ε) is cross-section; f(ε) is electron energy distribution function. Townsend coefficient of the ionisation reaction is where μ e is electron mobility. When the air consists of 80% N 2 and 20% O 2 , the Townsend coefficient of the ionisation reaction in air can be expressed as Attachment reaction mainly occurs between the electron and the O 2 molecule, and the corresponding Townsend coefficient can be obtained in the same way Townsend coefficient for the ionisation and attachment reaction is shown in Fig. 2.

Critical boundary condition
The influence of dielectric on discharge process is mainly considered in two aspects: (i) charged particles react on insulation surface and contribute secondary electrons; (ii) charge accumulates on the insulation surface, changing the distribution of the electric field. According to the treatment of surface reactions of insulating media in the literature [26,27], it is assumed that when the ions reach the wall, they transfer charge to the dielectric surface, charging the surface of the dielectric, and the charged particles transform into neutral particles.

Wall boundary condition:
The wall boundary condition is used when the plasma is in contact with a solid surface. At the wall, the exchange of electrons occurs through thermal motion, electro-migration and secondary electron emission. Compared with pure air gap discharge simulation, the surface discharge process needs to consider two aspects of secondary electron emission: (i) the secondary electron emission of the cathode caused by the positive ion striking the high voltage needle electrode; (ii) secondary electron emission caused by high-energy particles bombarding the material surface. The flux boundaries for electron density and electron energy density are as follows: The first term on the right-hand side of (16) and (17) refers to the changes in the number and energy of electrons caused by the disappearance/rebound of electrons on the wall due to thermal motion, where r e is the reflection coefficient. The value ranges from 0 to 1 and is set to 0 in this paper. v e,th is the thermal velocity; the second term is the electro-migration term, and the value of coefficient α e depends on the direction of electric field. Considering that the pulse electrical stress applied in this paper is negative polarity, the value α e of is 1 at the bar electrode and − 1 at the needle electrode; the third term is the electron quantity and energy gain caused by secondary electron emission, where γ is the second electron emission coefficient and ɛ is the average energy of second electrons. Based on pure air gap discharge simulation [28], the cathode secondary electron emission coefficient is set to 0.004. According to description of secondary electron from material surface reference DBD simulation [11], secondary electron emission of PI surface is set to 0.5; the average energy of secondary electrons ɛ = 2.5 eV [29]. The flux boundary for species k considers a flux contribution due to migration and an inward or outward mass flux determined by surface reactions. It is assumed that there are i = 1, 2, …, I kinds of surface reactions and k = 1, 2,…, K kinds of heavy species, then The rate constant of reaction i is expressed by the following equation: where Γ s is the total surface site concentration; γ i is sticking coefficient of reaction i, which is between 0 and 1; M k , c k , μ k and Z k are, respectively, the molar molecular mass, molar concentration, mobility and charge number of species k; sqrt(8RT/ πM k ) is the thermal velocity of species k; R is the universal gas constant; T is the gas temperature; σ k is potential characteristic length; v ki is the stoichiometric matrix; q i represents the surface reaction rate for reaction i and R surf,k is the surface rate expression for each species which comes from summing the surface reaction rates multiplied by their stoichiometric coefficients over all surface reactions.

Surface charge accumulation:
In surface discharge process, the surface of insulating medium will accumulate electric charge, resulting in discontinuity of electric displacement vectors on both sides of the dielectric interface, and the surface charge accumulation boundary condition [30] is D 1 and D 2 are, respectively, the electric displacement vectors on both sides (air side and medium side) of the insulation surface; σ s is the surface charge density, which can be computed by solving the following equation on the surface: where n·J i is the normal component of the total ion current density at the wall and n·J e is the normal component of the total electron current density at the wall.

Mesh generation and solution method
The model proposed above is solved by simulation software -COMSOL Multiphsics. The grid is divided into 146,543 cells. To ensure the computational speed and accuracy, the mesh should be non-uniformly divided. Near the air-PI interface, a very fine resolution is employed to resolve the steep gradients, as required. After many attempts, it is found that the space charge mainly concentrates in a thin layer 1 mm away from the dielectric interface in discharge processes, so mesh refinement is carried out in this area, and the cell size is controlled below 30 μm. The larger mesh size is used in the area which has less influence on the calculation results, and the maximum grid cell is not more than 150 μm. To avoid the negative value of particle density in the process of solution, the logarithmic finite-element method is used. The variables to be solved in the partial differential equation are logarithmised to get the logarithm of particle density, and then index it to get the actual particle density. Although the logarithmetics of variables increases the non-linearity of the equation, it greatly increases the stability of numerical solution. The environmental conditions are set at 20°C and 1.013 × 10 5 Pa. Initial electron density is 10 10 m −3 , initial positive and negative ions densities are both 10 17 m −3 .

Surface discharge under single pulse
A large number of simulation results show that the different repetition rates, peak values and pulse widths of the pulse applied by the needle electrode have some influence on the discharge propagation velocity. However, in general, there is no significant difference in the trend of surface discharge under repetitive pulses and the temporal and spatial evolution of electric field strength and electron density during discharge development from single pulse discharge with large amplitude and short rise time. Plasma channel always develops from needle to bar electrode gradually. Therefore, a single pulse with large amplitude and short rise time is selected to simulate the dynamic evolution of surface discharge, which is more intuitive and concise. The applied pulse voltage is shown in Fig. 3, with a peak value of 20 kV and a pulse width of 60 ns. The bar electrode is grounded and the rest of the boundaries are insulated.

Electric field distribution and particles dynamics behaviour
In the process of surface discharge, the distribution of electric field strength and electron density can visually show the development process of discharge, which is of great significance to understand discharge mechanism. It can be seen from Fig. 4 that the maximum electric field appears near the needle electrode in the initial discharge stage. Due to the strong electric field, a large number of electrons are ionised in the area near the needle electrode and the electrons move towards the bar electrode driven by electric field forces. As the speed of electrons is faster than that of positive ions, positive ions will stay in the ionisation region, so that the electric field between the needle electrodes and positive ions will be strengthened and collision ionisation near needle electrodes is intensified, thus, the development of discharge is promoted.
Figs. 5 and 6 show the dynamic distribution of the electron density and electric field strength along insulation surface. According to the effect of space charge on electric field distribution, the whole discharge process can be divided into three stages: (i) the development stage of electron avalanche (0-8 ns); in this stage, the effect of space charge on electric field can be neglected, and the number of electrons increases sharply from 10 10 to 10 17 m −3 of magnitude; (ii) the stage of electron avalanche converted into streamer (8-9 ns); in this stage, electron density growth has reached a plateau, but the influence of space charge on the external electric field is still small; (iii) streamer propagation (9-16 ns); in this stage, the space charge strongly distorts external electric field, and the area with the greatest electric field strength develops forward, which promotes the plasma channel on the insulation surface to extend gradually to the ground electrode, and finally penetrates the two electrodes at 16 ns.
It can be found that electron density and electric field have a good correspondence. With discharge propagation along the surface, the electron density and electric field strength show an increasing trend in the whole discharge area, and the maximum value of the electronic density always appears near the maximum of electric field strength, because the strong electric field leads to more intense ionisation and higher electron density. At any time, maximum value of the electron density and field strength appears near the needle electrode and the head of streamer, and increases gradually with the streamer propagation. What is the reason for the higher electron density near the needle electrode? On the one hand, the asymmetry of the needle-bar electrode results in the formation of a strong electric field near the needle electrode; on the other hand, the space charge further increases the electric field strength near the needle electrode. The composite field formed by the superposition of the external electric field and the space charge field reaches 10 kV/cm or so, and the strong electric field causes the intense collision ionisation near needle electrode. Similarly, there are two reasons for the high electron density at the head of the stream: (i) with the extension of the plasma channel, the insulation gap becomes shorter, and the applied voltage increases, thus the electric field strength increases; (ii) for the space charge distorts the external electric field, the electric field in front of the streamer head continues to increase, resulting in intense ionisation and increased electron density.
Distribution of positive and negative ions density in the development of surface discharge are shown in Fig. 7. As can be seen from the figure, the distribution law of positive ion density is similar to that of electron density, which is mainly because the main sources of both are collision ionisation. Compared with the positive ions, the negative ion density does not have a very obvious peak, but exists relatively uniformly in the plasma channel. Negative ions are mainly formed by the attachment of electrons and neutral molecules. Therefore, negative ions are abundant in areas with high electron density, but the negative ion density is generally 1-2 orders of magnitude smaller than that of positive ions. The reasons are as follows: (i) ionisation is the main reaction under strong electric field, resulting in a large number of positive ions; (ii) under the action of strong electric field, electrons move fast and are difficult to be trapped by oxygen molecules to form negative ions; (iii) the electric field strength in plasma channel is  Electron temperature is a characterisation of electron energy. The electron temperature distribution at different times in surface discharge is shown in Fig. 8. Electrons obtain energy from space electric field. The higher the electric field strength, the more energy is obtained by electrons, and the higher the electron temperature is. Therefore, the maximum temperature of electrons also appears in the head of the streamer. Electron temperature is roughly maintained at 5-6 eV in discharge process, and inside the plasma channel is lower, which is because electrons inside the plasma channel collide with other particles frequently, resulting in a large loss of electron energy.

Surface charge accumulation
The surface charge accumulation can change the spatial electric field distribution, which in turn affects the surface discharge process. As shown in Fig. 9, the insulation surface is concentrated with negative polarity charge. With the development of the discharge along PI surface, the amount of negative charge accumulation on the surface increases continuously, and the highdensity charge region gradually develops from the needle electrode to the ground electrode. The average charge density accumulated on the PI surface in the development process of surface discharge is shown in Fig. 10. It can be seen that the surface charge accumulation rate increases rapidly in streamer propagation stage (after 9 ns), and the average surface charge density reaches − 1.11 μC/m 2 at 16 ns.
To verify the reliability of the model, the air-PI surface discharge experiment was carried out under atmospheric pressure by using a needle-bar electrode structure. As the discharge process is very rapid under the simulation conditions, it is difficult to accurately capture the discharge process in the existing experimental conditions. The surface charge distribution is difficult to realise real-time online measurement. In order to accurately capture the discharge process and obtain the change of the surface charge distribution with the development of the discharge, the adjustment voltage is selected in the experiment to slow down the discharge development process. The whole period from the start of discharge to the flashover is defined as the surface life, and it is divided into the first to fourth discharge development stages by the 20, 60, 80 and 100% nodes of surface life, corresponding to 3.2, 9.6, 12.8 and 16 ns moments in simulation. The development of surface discharge was recorded by a camera, and surface charge accumulated in different stages of surface discharge were measured off-line by using a feedback electrostatic probe and a Trek Model 347 electrostatic voltmeter. As can be seen from Fig. 11, as the continuous development of surface discharge, the plasma channel gradually extends from the tip of the needle electrode towards the bar electrode along the surface of PI film, and finally penetrates the two electrodes to form a flashover arc. Also, the discharge trajectory is approximately linear, which verifies the rationality of the dimensionality reduction of the model. Fig. 12 shows the PI surface charge accumulation at different stages of the surface discharge (four stages shown in Fig. 11). The surface charge measurement range shown in Fig. 12 is a square area with 1.5 cm of side length near the needle tip. It can be seen that the PI surface area has a negative polarity charge and there is a high density  In general, the terminal current is due to the conduction current in a dielectric medium, and the displacement current is due to the time rate of change of the surface charge on the electrodes, as expressed by the following equation, which is known as the integral Ohm's law: As shown in Fig. 13, the calculated current shape is close to the measured current pulse, and the calculated current peak value is slightly 4% lower than the measured value.

Surface discharge propagation velocity
Surface discharge propagation velocity is closely related to the electric field strength in front of plasma channel. The stronger the electric field, the faster the discharge propagation velocity [31]. The electric field strength in front of plasma channel is formed by the superposition of the background field strength and the additional field strength produced by space charge and surface charge. Fig. 14 shows tangential electric field strength in front of plasma channel at different time. With the development of surface discharge, the tangential electric field strength and the growth rate are increasing. Taking arbitrary time t 0 and recording the position of the streamer head as s 1 and s 2 , respectively, at t 0 −Δt and t 0 + Δt, the propagation velocity v = (s 2 −s 1 )/(2Δt). Observing the variation of the velocity with time, it is found that in the development stage of the electron avalanche, the discharge propagation velocity along the surface is smaller than 10 5 m/s; in the stage of electron avalanche converted into streamer, the velocity increases more rapidly, the magnitude is 10 5 m/s; in the stage of streamer propagation, the velocity of the early stage (9-14 ns) increases steadily, and then increases sharply after 14 ns. It reaches 10 6 m/s, which is consistent with the simulation results in [32]. Fig. 15 shows more intuitively the relationship between the surface discharge propagation velocity v and the tangential electric field E t in front of plasma channel. The fitting function used is v = 2.975 × 10 −16 × E t 3.518 . The experimental results in [33] show that the component along the surface of streamer velocity increases significantly with the increase of applied electric field, which is consistent with the simulation results in this paper.

Inflection point caused by pulse repetition rate
The peak value of pulse voltage is fixed as − 1.2 kV. Figs. 16a-f show the electron density distribution when the plasma channel penetrates the two electrodes under pulse voltages with different repetition rates. It can be seen from Fig. 16 that the maximum electron density always appears in the streamer head, and the electron density ranges from 10 16 -10 17 m −3 . In addition, the simulation results show that the relationship between the discharge propagation time and the pulse repetition rate is not monotonous, as shown in Fig. 17, they are of 'u' curve relationship: with the increase of repetition rate, it first decreases and then increases, and there is a inflection point near 1667 kHz. This may be due to the fact that when the repetition rate of pulse voltage is very high, the voltage changes rapidly, creeping process cannot keep up with the change of voltage, and the plasma channel cannot penetrate the two electrodes in one cycle. The electron turn-back phenomenon appears in the discharge process with the change of pulse voltage, which makes the formation time of plasma channel longer. When the discharge velocity is fast enough that plasma channel can penetrate two electrodes in one cycle, the discharge propagation time is mainly determined by the thermal effect. The lower the repetition rate, the longer the time required for the formation of plasma channels through the two polarities. To confirm above explanation, two repetition rates (f 1 = 1250 kHz, and f 2 = 5 MHz) in front and back of inflection point were selected to observe the surface discharge process. As shown in Fig. 18a, when the pulse repetition rate is lower than that of inflection point, the plasma channel continuously develops from the needle electrode to the ground electrode. It can be seen from Fig. 18b that electrons turn back between 300 and 400 ns, and the maximum electron density moves gradually from the head of the streamer to the vicinity of the needle electrode, which is because the lower applied voltage at this time resulting in smaller external electric field E a , while the previous surface discharge produces a large number of positive and negative charge space separation, which results in the formation of space charge field E sc in the discharge channel (see Fig. 19). The internal E sc of the streamer channel is opposite to the direction of E a , which weakens the external electric field. When the E a decreases continuously, the compound field strength formed by E a and E sc in the plasma channel is opposite to the applied electric field E a , the electrons move to needle electrode driven by compound field, the number of electrons in the front of plasma channel decreases. Due to the weak external electric field, the electric field in front of plasma channel is insufficient to maintain the further propagation of the surface discharge. Until the applied voltage increases, the composite electric field strength can reach the required field strength for the development of surface discharge and the surface discharge can further advance to the ground electrode. After several cycles, a plasma channel through the two electrodes is formed and the discharge process is completed.

Influence of pulse voltage peak and pulse width on discharge development process
The fixed pulse repetition rate is 10 MHz, and changed voltage peak V a = −1.2, − 1.6, − 2, − 2.4, − 2.8 and − 3.2 kV. It is found that the relationship between the peak value of pulse voltage and the discharge propagation time is shown in Fig. 20. The larger the peak value of voltage, the faster the discharge will develop, which is because the higher voltage, the stronger the electric field, the more intense the discharge, the more electrons, the easier the plasma channel will develop; and as the increase of the voltage peak value, the propagation time decreases at a slower and slower rate and eventually becomes saturated. Fig. 21 shows the effect of pulse width on the discharge process. When the peak value and repetition rate of the pulse remain unchanged, the larger the pulse width is, the faster the discharge develops and the shorter the time required for plasma to penetrate the two poles is. The effect of pulse repetition rate on discharge propagation time is shown in Fig. 22 at different peak values of pulse voltage. It can be seen from the figure that at any repetition rate, with the increase of peak voltage, the discharge development accelerates and the plasma penetration time shortens, which is consistent with Fig. 20. At any voltage, there is 'u' curve relationship between the pulse repetition rate and the discharge propagation time, which is consistent with the conclusion in Section 4.1. Further analysis shows that with the increase of peak voltage, the pulse repetition rate corresponding to the inflection point increases gradually. From Section 4.1, it can be known that the occurrence of inflection point is related to the discharge propagation velocity. When the pulse repetition rate increases to the point where the discharge channel cannot be completed in one pulse cycle, the inflection point will appear. The increase of applied voltage will lead to the acceleration of plasma channel development, and the time of penetrating the two electrodes is shortened, thus the pulse period at inflection point decreases and the repetition rate increases correspondingly.

Conclusion
In this paper, needle-bar electrode structure is used to simulate surface discharge process, and a non-equilibrium plasma model of air-PI surface discharge is established. On the basis, the dynamic development process of surface discharge under pulsed electric stress is implemented, and spatio-temporal evolution of electric field, particle density, space and surface charge density and other microscopic parameters are obtained. Based on the simulation under the needle-to-bar electrode with 1 cm gap, the following conclusions can be drawn: (i) As air is a mixed gas, it is extremely difficult to establish a complete plasma reaction system. A simplified set of reactions is used to describe chemical reactions in the process of air discharge, which greatly reduces the difficulty of modelling, and makes it possible to simulate the discharge along the surface with a 1 cm gap.
(ii) The surface discharge process under single pulse can be divided into three stages: electron avalanche, electron avalanche converted into streamer and streamer propagation. At the initial stage of discharge, the number of electrons increases rapidly and concentrates near the needle electrode. With the development of discharge, a large number of electrons are produced at the streamer head, and the plasma channel extends to the ground electrode. Finally, when the two electrodes are penetrated, the electron density reaches magnitude of 10 20 m −3 . (iii) The development trend of surface discharge under repetitive pulses is basically the same as that under single pulse: The plasma channel extends from needle electrodes to ground electrodes continuously. Discharge propagation time and pulse repetition rate are of 'u' curve relationship, and the inflection point gradually moves towards the high repetition rate direction with the increase of the peak pulse voltage. (iv) The numerical model established in this paper can be used to analyse the discharge process of centimetre level of gap along airsolid insulation surface under nanosecond pulses, which can make up for the shortcomings of experimental methods. However, how to apply the numerical simulation methods in larger time scale (μs level or even ms level) and spatial scale (decimetre level or even metre level) is worthy of further study.

Acknowledgment
This work was supported by National Natural Science Foundation of China (grant no. 51929701).