Automatic Optimization of Gas Insulated Components Based on the Streamer Inception Criterion

: Gas insulated transmission lines (GILs) are used in electrical systems mainly for power transmission and High Voltage substation interconnection. In this paper, we focus on the development of complex numerical tools for the optimization of gas insulated HVDC components by the estimation of realistic electric ﬁeld distribution and the voltage holding of the designed geometry. In particular, the paper aims at describing the correct modelling approach suitable to study high voltage components in DC, considering the nonlinear behaviour characterizing the electrical conductivity of solid and gas insulators. The simulated ﬁeld distribution is then adopted to estimate the voltage holding of the dielectric gas, with a convenient engineering technique, based on the streamer criterion. These two tools are integrated in an automatic optimization package developed in COMSOL ® and MATLAB ® , with the purpose of adjusting the critical geometry features, suffering from excessive electrical stress and possibly giving rise to electrical breakdown, in order to guide the designer towards a robust solution.


Introduction
High Voltage Gas Insulates Lines (GILs), because of their compact size, large-capacity transmission and stable operation, have become of crucial importance in many applications. Currently, most of the applications are AC GILs, but the improved performance of power electronics used in conversion stations and the increasing demand for electricity transport in some cases for hundreds of kilometres because of the energy turnaround, HVDC GILs are increasingly attracting attention [1][2][3]. Usually, SF 6 at 0.6 MPa is used as insulating gas despite its high Global Warming Potential (GWP), about 23,500 higher than CO 2 . In the last years, solutions with lower impact on the environment, like SF 6 mixtures and fluoronitriles were taken into account [4][5][6].
The conceptual design of GIL must be defined aiming at reducing their complexity and economic impact. Conflicting demands must be met: reduce gas pressure, minimize mechanical stress and gas volume, minimize geometrical dimensions (e.g., cross section), in order to reduce cost, weight and impact on the route footprint or building design, by guaranteeing at the same time the dielectric strength of the system avoiding breakdown events. Moreover, the main issue affecting HVDC-GIS systems, and not HVAC-GIS, is related to electric charge accumulation on the interface between gas and solid insulators, during the charging process from the initial capacitive field to the stationary resistive distribution. Charge accumulation on dielectric surfaces may lead to a decreasing flashover voltage [7]. Such decrease is more pronounced for DC fields than AC fields [8]. However, different factors affect the flashover voltage under AC and DC conditions [9].
Firstly, the present work describes the numerical model used to correctly simulate the electric field distribution in HVDC-GIS systems, self-consistently taking into account

Electric Field Modelling
The electric field distribution within insulating materials under DC voltage is different from that under AC excitation, depending both on the electrical permittivity and the conductivity of materials [10,11]. In order to correctly study the charge accumulation mechanism in HVDC-GIS, the non linear conductivity in the solid dielectrics, the physical processes affecting the motion of charge particles in the gas and the surface conductivity at the gas-solid interfaces must be considered at the same time [12].

Gas Domain
Differently from solids, where it is more appropriate to speak about the concept of electrical conductivity, which is usually a space function depending on the temperature and possibly on the electric field itself, a gas model is required for the correct analysis of electric fields under DC voltages. The model takes into account all the physical mechanism affecting the electrical conduction in the gas: generation, recombination, motion and diffusion of charge carriers [13]. In electronegative gases like SF 6 , a negligible amount of free electrons is present in normal conditions, and the motion of positive and negative ions is governed by a set of convection-diffusion Partial Differential Equations (PDE) [14]: where n + and n − are the ion number densities, ∂n IP /∂t is the ion-pair generation rate, R is the recombination coefficient and µ + , µ − the ion mobilities, and D ± is the diffusion coefficient related to mobility via Einstein relation [15]: with k B , T and e the Boltzmann constant, temperature, and elementary electric charge respectively. Equations (1) and (2) must be coupled with the Poisson equation to calculate a self-consistent electric potential φ G in the gas domain: while Gauss law allows to impose the ion densities as the source of the electric field in the gas domain ∇ · (ε G E G ) = e(n + − n − ), with ε G = ε 0 ε r,G the permittivity of the gas. Thus, the ion drift due to the applied electric field and the diffusion due to local differences in the charge density [16], contribute to the gas current density, resulting into the following: Equations (1)- (5) are combined in the implemented gas model characterized by a strong nonlinear behaviour.

Solid Domain
In the solid (bulk) insulator domain, in the general case of transient phenomena, the Gauss' law, the equation of continuity, and Ohm's law hold: where k I and ε I are the electrical conductivity and the permettivity of solid insulator and ρ is the bulk charge density. The transient equation for the space charge density in solid insulator is derived combining Equations (7)-(9): In deriving Equation (10), the bulk permittivity is assumed to be constant, since it is only slightly affected by temperature and electric field intensity. On the contrary, the solid insulator electrical conductivity is generally considerably affected by electric field and temperature distributions, i.e., k I = k I (E, T). The following relation can be considered [17] for epoxy insulators: with a = 0.08 mm/kV, b = 0.1°C −1 and 10 −20 < k I,0 < 10 −14 S/m. For the subsequent analyses the temperature was assumed to be uniform and the solid conductivity weakly dependent on the electric field, therefore a constant value is adopted since this doesn't affect the main outcomes presented.

Gas-Solid Interface
In addition to the bulk current flowing in the solid and gas insulators, also the current flow due to surface conductivity k S and the tangential component of the electric field E t along the insulator surface contributes to the process of accumulation of surface charge density σ S at the gas-solid interface, during the capacitive-resistive transition. The charging process is described by the thin layer continuity equation [10,14]: where J I,n and J G,n are the normal components of the current densities in the solid insulator and in the gas, respectively.

Boundary and Initial Conditions
The solution of nonlinear time dependent PDEs requires appropriate boundary and initial conditions. For the high voltage and low voltage electrodes, Dirichlet boundary conditions are specified for the electric potential. On gas-solid interfaces, the continuity of electric potential is imposed: For the ion advection-diffusion Equations (1) and (2), denoting with n G the outward normal vector of gas, the following boundary conditions are assigned [18]: • No positive ions are emitted from surfaces where n G · E G < 0, i.e., n + = 0 and no diffusion of negative ions occurs from these surfaces, i.e., −D − ∇n − = 0. • No negative ions are emitted from surfaces where n G · E G > 0, i.e., n − = 0 and no diffusion of positive ions occurs from these surfaces, i.e., −D + ∇n + = 0.
The time dependent problem also needs initial conditions. In this case, zero electric potential and equilibrium state for the ion number densities are assumed for the whole system [10]: In order to achieve accurate results, the level of mesh refinement in specific areas and features such as triple points and sharp edges are critical. Particularly in case of time dependent 3D analyses, a uniform refinement all along the domain is wasteful in terms of computation time and resources. Local refinement is of course convenient and a so-called goal-oriented adaptivity based on a quantity of interest which is function of the solution, as described in [19], is of particular interest in this case. Indeed, as will be described in Section 4, the objective function adopted for the optimization algorithm might be used in an inner loop of the optimization algorithm in order to improve the solution of the electric field distribution.

Streamer Breakdown
Knowledge of dielectric breakdown of a gas insulated system requires modelling the basic precursors, such as the streamer inception, streamer propagation, and the streamer to leader transition. With the latter appearing when the streamer discharge reaches a certain length [20]. From the engineering point of view, it is important to highlight geometrical features suffering from excessive electrical stress in order to optimize their shape. Knowledge of the electric field distribution alone, thus considering the areas where the field exceed a critical value, can be misleading. As a matter of fact, the breakdown in order to take place requires the fulfilment of particular conditions along the whole distance from cathode to anode.
The streamer inception criterion for electronegative gases, e.g., air, N 2 or SF 6 is based on avalanche theory [21]. When an electron avalanche reaches a critical size N cr in the streamer head, the streamers can result into a partial or complete breakdown of insulating gas. The typical value of N cr is about 10 8 electrons [22]. The streamer breakdown criterion reads [23]: whereᾱ(E, P) [1/m] is the electric field strength E and gas pressure P dependent effective ionization coefficient [24]. It is defined as the difference between ionization and attachment coefficients, i.e.,ᾱ := α − η, both of which are function of the electric field. The quantityᾱ/P = f (E/P) is often used in this context and is called reduced effective ionization coefficient. The integration (15) is performed along electric field lines Γ whereᾱ > 0, from the point of maximum field and ending whereᾱ ≤ 0. To take into account this constraint, the integrand in (15) is multiplied by the Heaviside function Θ(ᾱ) [25].
When dealing with electrode-insulator systems, the flashover voltage, corresponding to occurrence of the electric discharge between the electrodes, can be estimated using empirical relations or analytical formulas derived from a circuit model [20].
In this paper, we focus only on the first precursor of breakdown phenomena, i.e., the streamer inception. Although this criterion does not provide actual information on the discharge and flashover, it allows us to obtain information on the critical areas of the geometry in a simple and relatively low computational way.

Effective Ionization
A key point of this work is the calculation of the effective ionization. As a matter of fact, the computation of streamer integrals (15) requires knowledge ofᾱ for a large range of operational parameters. For common insulating gases empirical relations well-fitting experimental data over a wide range of E/P are listed in the literature [21,26]. Using the ideal gas law, relating gas pressure and density N [1/m 3 ], the reduced effective ionization coefficient of SF 6 , can be expressed as [26]: for 360 < E/N < 5000 townsend (Td). For given electric field values, an increase in gas pressure reflects in a decreasing trend ofᾱ, thus reducing the probability of streamer inception. As stated in the introductory part, the use of alternative insulating gases, such as SF 6 /N 2 mixtures [27] or fluoronitriles, requires knowledge of voltage holding capabilities [28]. Using the streamer inception criterion to estimate such voltage holding capability, requires an expression for the effective ionization coefficient.ᾱ can be obtained experimentally or, more conveniently, can be calculated numerically solving the collisional Boltzmann equation for a specified set of electron collisional reactions [29,30]. This procedure starts with the definition of the set of reactions, R i , each characterized by its scattering cross section σ i ( ), function of the incident electron energy. The electron energy distribution function (EEDF) f ( ), according to the kinetic theory of plasmas, satisfies the Boltzmann equation, which describes its evolution in a six-dimensional phase space [31]: where C[ f ] is the collisional term. From the electron energy distribution function, the rate coefficient for the electron collision reaction R i can be calculated [32]: with γ = √ 2e/m e . The swarm coefficient α i /N, corresponding to the collision reaction with cross section σ i ( ), is directly computed using the relationship [32]: where W is the electron drift velocity.

Shape Optimization
In this section, on the basis of the previous considerations, the shape optimization procedure is presented. It is constituted by three main tools: the electric field solver, the streamer breakdown detector and the numerical optimization algorithm. The problem of minimizing a given objective function (OF) f (x) under D design variables described by the vector x, can be formally defined as: To solve such optimization problems (20), different techniques can be considered [33,34]. Here the heuristic Differential Evolution (DE) algorithm is adopted, a parallel direct search method presented by Prince and Storn [35], whose scheme, shown in Figure 1, is based on typical evolutionary schemes involving mutation, crossing over and selection. The DE scheme is used in several industrial applications, e.g., for the optimization of Near Field Communication (NFC) antennas [36][37][38]. The basic principles of DE algorithm are hereafter summarized. The initialization procedure creates a random population of N individuals in a D-dimensional parameter space, within given boundary constraints. Each generation k deals with: parameter vectors. The mutation operation expands the search space adding difference vectors (donors) to a base individual: where v i i ∈ {1, . . . , N} is called the donor vector, while first x on the right-hand side is the target. F represents the mutation probability factor. The indexes r 1 , r 2 , r 3 ∈ {1, . . . , N} have to be chosen mutually different and also different from the index i [35]. Usually, best performances are obtained choosing x r1 = x best [39]. In the recombination phase, parameters from the target and donor individuals are mixed in order to form the so-called trial individuals: where CR is the crossover ratio, rand i,j random numbers and irand a random integer. Vectors u i = {u i,j } j = 1, . . . , D, are called trial vectors. The last block of Figure 1 is the selection. Here, a greedy scheme is applied in order to select the next generation: In this work, the automatic shape optimization proceeds as follows. The area of interest is parametrized using, e.g., a Bèzier curve or spline for which its control points represent the design variables and an initial mesh of the domain is generated. In this work, the (stationary) electric field and the voltage holding predictor are then computed according to the previously described methods; the objective function is the maximum value of the effective ionization coefficientᾱ. The objective function is therefore evaluated and new design variables are automatically generated by the search algorithm. The procedure repeats until the termination criterion is reached. Here, the stop criterion is prescribed by the maximum number of search algorithm iterations it max .
When the number of objective functions to be minimized in (20) is greater than one, the optimization problem becomes a multiobjective optimization problem. Such problems can be solved using again evolutionary algorithms (DE), based on Pareto optimal solutions [39]. For these problems, the aim is to find compromise solutions instead then a single optimal one as in scalar optimization problems, e.g., using the method of objective weighting [17].

Electric Field Modeling in Air
The model described in Section 1 has been applied to determine the HVDC-GIS electric field distribution for the geometry shown in Figure 2. It consists in a 2D axisymmetric solid insulator with radius r I = 27 mm, height h I = 150 mm with embedded metallized cavities. The solid insulator separates a High Voltage (HV) electrode from a Low Voltage (LV) one at the bottom, all immersed in gas with a certain pressure in a cylindrical chamber with height H = 600 mm and radius W = 300 mm. The two feedthroughs are electrically insulated from the chamber at ground potential through HV and LV bushings. In this test-case, air at atmospheric pressure was considered as dielectric gas.  The voltage applied to the high voltage electrode is U DC = 1 kV, while the lower voltage is referred to ground. Epoxy resin with ε r,I = 5 and k I = 3.33 × 10 −19 S/m is chosen for the solid insulator. Following [10], for atmospheric air, the ion mobilities are set to µ + = 1.36 × 10 −4 m 2 /(Vs), µ − = 1.87 × 10 −4 m 2 /(Vs), the recombination coefficient R = 1.4 × 10 −12 m 3 /s and the ion-pair generation rate ∂n IP /∂t = 7 × 10 6 IP/(m 3 s). In this example, the surface conductivity of gas-solid interface is set to zero, i.e., k S = 0. The electric field is numerically solved with the commercial FEM software COMSOL ® Multiphysics. Equations from (1) to (5) are implemented in the gas domain, Equations (10) and (11) in the solid domain. The surface charge temporal variation described by Equation (12) is implemented separately. The condition for calculating the stationary-resistive state are: and The simulated electric potential and charge density distributions along the creepage distance on the gas-solid interface, at different time instants, are reported in Figures 3 and 4. Figure 5 shows the temporal variation of the electric field components normal and tangential to the insulator surface.

Streamer Inception in SF 6
The numerical implementation of voltage holding predictor described in Section 3 is hereafter applied to pure SF 6 , for which a simplified set of electron-molecule reactions is listed in Table 1 [30]. The collisional cross section data σ i ( ) are taken from LXCAT online database [40]. (17) is computed in the Plasma module of COMSOL ® Multiphysics. The reduced Townsend coefficient α/N and the reduced attachment coefficient η/N are automatically obtained in the post processing analysis. Figure 6 shows a comparison ofᾱ/N calculated with the aforementioned method and the empirical relation (16). The results are in good agreement.
The geometrical model of Figure 7 is used to determine the voltage holding of SF 6 at 0.1 MPa, with an applied voltage of U DC = 250 kV. The proposed electrode-insulator system differs from GIL configurations, but corresponds to typical experimental arrangements adopted for investigations on GIS applications. For the electric field calculation, stationaryresistive conditions are taken. From the literature, it is known that SF 6 ion mobilities depend, on pressure as [41]: therefore, in this example, the values of µ ± = 3.6 × 10 −5 m 2 /(Vs) are taken. The recombination coefficient and ion-pair generation rate are R = 1.7 × 10 −12 m 3 /s and ∂n IP /∂t = 1 × 10 7 IP/(m 3 s), respectively [42].
Once the electric field is computed, streamer criterion can be verified. The effective ionization coefficient previously computed is interpolated to obtain its value for specific electric field values from the simulation results.
Streamer integrals are calculated with COMSOL ® Charge Particle Tracing (CPT) module, integrating (15) along particle trajectories and using a massless formulation [43,44]. The results are shown in Figure 7. It can be observed that, for the assigned operational parameters, the streamer breakdown inception occurs: there exist field lines, originating from the high voltage electrode, for which K > K cr . Moreover, the peak value of the simulated electric field is above 8.9 kV/mm, corresponding to the breakdown electric field at 0.1 MPa [45].  (15) satisfying streamer criterion for discharges originating from high voltage electrode. The peak value of simulated electric field is above the breakdown electric field of SF 6 of 8.9 kV/mm at 0.1 MPa. The inset detail shows the region whereᾱ > 0 in proximity of the high voltage electrode.

Shape Optimization of High Voltage Electrode
In this test-case the objective function to be minimized is the maximum value of the effective ionization coefficient; in particular we are searching for a new shape of the high voltage electrode, for which the streamer criterion (15) is not satisfied. The optimization tool described in Section 4 has been developed using COMSOL ® with MATLAB ® LiveLink TM . The differential evolution algorithm was written in MATLAB ® and its main parameters have been set as follows:  Figure 9 shows the flowchart of the optimization, highlighting the steps executed in COMSOL ® and MATLAB ® , while Figure 10 displays the electric field map for the optimized configuration. Starting from the peak value of 10.5 kV/mm, the optimization procedure reduces the maximum electric field to about 8.5 kV/mm, below the critical value for streamer inception, as it can be seen in Figure 10 where all calculated trajectories are associated to a zero value for the streamer integrals (15).    (15) and electric field map for for the optimized electrode shape. The peak value of the electric field is reduced to about 8.5 kV/mm and the streamer criterion is not exceeded for any of the calculated trajectories for which the integral (15) is always zero.

Shape Optimization of Insulator
The optimization tool described in Section 4 can be also used to mitigate the surface charge accumulation on dielectric surfaces, that, in industrial applications, is determined using the Charge Inversion Algorithm [46,47], which can be accelerated by exploiting hierarchical matrices as shown by many works in the literature [48][49][50]. As reported in [7], the surface charging depends on normal and tangential components of the electric field on the solid insulator at the initial capacitive state [51]. With the purpose of reducing the charge accumulation and mitigate the probability of surface flashover, the field components can be minimized through a shape optimization of the insulator [52]. In this case, the minimization problem can be defined as a multiobjective optimization problem [53], where the design variables define the shape of the insulator surface and the objective function becomes a vector constituted by the normal E n and tangential E t (or its gradient ∇E t ) components of electric field on solid insulator, i.e., The multi-objective optimization has been applied to the layout of the insulator showed in Figure 11 (dashed blue). The shape of insulator surface is described with a Bézier curve with seven control points, of which only P3 and P4 are used as design variables, due to the symmetry of the system, i.e., the parameter vector became x = [r P3 , r P4 , z P3 ]. P2 (and P6) are kept fixed in order to obtain reasonable shapes. The cross-section profile for the optimized geometry is presented (continuous red) together with the profile of the original configuration. The related electric field components E n , E t of the two configurations are shown in Figure 12.

Conclusions
A numerical tool, implemented in COMSOL ® and MATLAB ® , for the automatic optimization of HVDC-GIS, combining the auto-consistent modelling of electric field distribution in solid and gas dielectrics and exploiting streamer inception criterion to evaluate the voltage holding of a predefined geometry, has been presented. The tool was firstly applied for the transient analyses of an air-insulated gas-solid system and subsequently for the shape optimization of the high voltage electrode in a pressurized SF 6 chamber. Moreover the tool allows the possibility of optimizing the shape of the insulator itself with the purpose of minimizing the electric field components, in order to reduce the effect of accumulated charge, which can drastically reduce the flashover voltage. The automatic optimization of the critical details of the designed geometry is a crucial aspect from the engineering point of view. The possibility to robustly estimate the effective ionization coefficient of different gases, such as pure gas and gas mixtures for which there are no experimental data in terms of effective ionization coefficient, and in different pressure conditions is a key feature of this tool, which allows a widely application of the proposed method.