Investigating the E ﬀ ect of Several Model Conﬁgurations on the Transient Response of Gas-Insulated Substation during Fault Events Using an Electromagnetic Field Theory Approach

: Assessment of very fast transient overvoltage (VFTO) requires good knowledge of the behavior of gas-insulated substation when subjected to very high frequencies. The international standards and guidelines generically present only recommendations regarding the VFTO suppression without a technical and mathematical background. Therefore, to provide an accurate image regarding the critical locations across a gas-insulated substation (GIS) from a transient response point of view, a suitable modeling technique has to be identiﬁed and developed for the substation. The paper aimed to provide an accurate assessment of the GIS holistic transient response through an electromagnetic ﬁeld theory (EMF) approach. This modeling technique has always been a di ﬃ cult task when it came to gas-insulated substations. However, recent studies have shown that through suitable Computer-aided design models, representing the GIS metallic ensemble, accurate results can be obtained. The paper investigated several simpliﬁcations of the computational domain considering di ﬀ erent gas-insulated substation conﬁgurations in order to identify a suitable modeling approach without any unnecessary computational e ﬀ ort. The analysis was performed by adopting the partial equivalent element circuit (PEEC) approach embedded into XGSLab software package. Obtained results could provide useful hints for grounding grid designers regarding the proper development and implementation of transient ground potential rise (TGPR) mitigation techniques across a gas-insulated substation.


Introduction
Typically, the gas-insulated substation (GIS) requires 10-25% of the area allocated to a conventional air-insulated substation due to its superior dielectric characteristics of the insulation material, sulphur hexafluoride (SF6), compared to air. SF6 gas is widely used in high-voltage energy applications also because of its high electric arc quenching capacity, which is chemically inert, nontoxic, and nonflammable [1]. Due to dielectric material strength reasons, local values of the electric field inside GIS enclosure, there are two types of GIS topology: three-phase system encapsulated in a single coaxial pipe for nominal voltages below 145 kV and each phase conductor individually encapsulated in a coaxial pipe for nominal voltages above 145 kV [2]. Besides the occurrence of the widely known and studied short circuit at power frequency fault, gas-insulated substations are subjected to severe transient regimes generated by the switching events. As a result of the relatively low speed of the disconnector switch contacts' movement during switching events, 0.6 cm/s [3], there is dielectric breakdown phenomena occurrence followed by the appearance of electric arc in the contact cavity causing the generation of transient overvoltage characterized by very-high-frequency (Very Fast Transient Overvoltage, VFTO) [4].
The amplitude and waveform of the very fast transient overvoltage (VFTO) depends on the configuration of the GIS substation, the trapped charge stored by the GIS equipment [5], the resistance of the electric arc formed between contacts [6], the capacity at the transformer terminals [7], and the gas pressure [7]. Measurements and digital simulations' amplitudes of transient overvoltage were determined to be between 0.5 p.u. and 2.05 p.u. [8][9][10]. During switching events in GIS, the transient electromagnetic wave travels toward the metal enclosure through sulphur hexafluoride (SF6)-air bushings and the metallic flanges located between each particular metallic enclosure [11,12]. Moreover, when the difference of the potential between the inner wall of the enclosure and the phase conductor exceeds the breakdown voltage of the dielectric material an electric arc is initiated, generating a short circuit between both metallic structures. When the voltage increases over the breakdown voltage of the insulating arrangement, an arc discharge takes place [13]. This is characterized by a heavy flow of current through the gas between the electrodes and the high dissipation of energy in the form of heat. The breakdown voltage of pure SF6 gas (below 1 bar gas pressure) is around 240 kV when a 20-mm gap between experimental electrodes is considered [14].
When a short circuit between the phase conductor and the inner wall of the enclosure occurs, due to voltage breakdown phenomena, an important amount of energy is leaked toward the metallic enclosure. Usually, GIS module enclosure is connected from place to place through grounding leads to the earthing system of the substation, thereby generating dangerous levels of transient ground potential rise (TGPR). In order to accurately compute and assess the TGPR occurring in GIS, the three-dimensional configuration of the metallic shell as well as grounding grid conductors need to be considered.
There are several locally applied mitigation techniques (like ferrite rings, shunt resistors, etc.) proposed in the literature [15][16][17]. However, these cannot be generally adopted due to mechanical constraints and they usually require high additional financial efforts. Moreover, the international standards and guidelines provide recommendations regarding the VFTO suppression without a technical and mathematical background [1]. An accurate assessment of the holistic transient response of the substation can highlight the critical locations across the GIS and its vicinity. As a result, suitable mitigation techniques can be applied only where they are necessary.
There are two main approaches widely used when GIS analysis needs to be employed: circuit theory and electromagnetic field theory approaches. The former is the predominant method adopted and used in industry and literature for simulation and investigation of the transient regimes in GIS. The circuit theory approach is based on the representation of encapsulated ensemble by equivalent electrical circuits and distributed parameters: propagation speed, equivalent capacitance, equivalent inductance, and characteristic impedance providing solutions in the frequency and time domain [18][19][20][21][22][23][24]. The circuit formulation methodology is based on Kirchhoff's laws embedded in different software interfaces providing numerical solutions for an energy system operating at a known voltage level, based only on the transverse mode of propagation of the electromagnetic wave, TEM, meaning that the method neglects the electric and magnetic field in the direction of propagation. [25]. The method neglects the effects of propagation losses (skin effect, etc.), which result in a lower damping coefficient for very-high-frequency components associated with transient regimes. An important drawback of the method consists of the impossibility to consider the three-dimensional configuration of the GIS metallic enclosure, as well as of the grounding grid conductors during the computational process. Considering such complex metallic structures possessed by a gas-insulated substation configuration (e.g., enclosure, structure resistance, grounding grid conductors) located in a relatively small air Energies 2020, 13, 6231 3 of 19 volume, electromagnetic couplings will occur between different metallic subsystems, which cannot be quantified in the final solution provided by circuit theory approach.
On the other hand, the electromagnetic field theory approach can be applied toward GIS three-dimensional geometries regardless of the complexity of the model to be studied [25]. Full-wave numerical electromagnetic analysis (NEA) methods are defined as methodologies allowing the direct numerical solution of Maxwell's equations in both frequency and time domain [25]. These methods are becoming the most promising approaches to study complex transient phenomena that cannot be straightforwardly solved by means of circuit theory or transmission line approaches (e.g., by using electromagnetic transient program (EMTP)). Computational electromagnetics are getting increasing attention not only in the research fields but also in the industrial fields as well as in gas-insulated substation analysis. A method solving Maxwell's equation directly can be classified into a differential equation-(DE) and an integral equation-(IE) based method [25]. The most commonly used numerical approaches in solving Maxwell equations, applicable to gas-insulated substation analysis, are method of moments (MoM) finite element method (FEM), finite differences time domain approach (FDTD) and partial equivalent element circuit (PEEC) method, which will be further applied and discussed. In the vast majority of the case studies available in the literature, MoM, the IE-based method approach applied to gas-insulated substation analysis, is focused on safety assessment during steady state, unbalanced condition, and power frequency faults [26][27][28]. Although the mentioned studies provide a comprehensive safety assessment analysis, the proposed Computer-aided design based (CAD) models representing the GIS metallic enclosure and adjacent inner phase conductors are not properly validated.
Current studies present FEM numerical approach as a viable solution for dielectric design, regarding GIS insulation material, in order to ensure that the electric field and its gradient are within acceptable limits, below critical values imposed by the manufacturer [2]. Furthermore, studies containing partial discharge detection in dielectric media (dielectric strength analysis) alongside mechanical stress analysis and voltage breakdown between disconnecting switch (DS) contacts are available in public literature [29][30][31][32][33]. FEM software packages can handle geometries regardless of the required accuracy (DS contacts' chamber). However, the limitation of the method arises when large geometries need to be included in the computational domain. FDTD, similar to FEM, is based on the DE form of Maxwell equations, resulting in the necessity to consider the entire volume space containing the geometry under study as computational domain, which implies a high amount of digital resources [34][35][36][37][38][39]. By adopting suitable boundary conditions (perfect matching layer [36], perfect electrical conductor [36]), the area of interest can be extended, considering several simplifying assumptions. Nevertheless, the computational effort required by considering the entire GIS metallic enclosure during the analysis represents an important drawback of the FDTD approach. The PEEC method is based on an integral equation description of the geometry that is interpreted in terms of circuit elements: partial inductances and partial capacitances, a partial potential, which represents the intermediate elements between electromagnetic field approach and interpretation of Maxwell's equations in circuit domain [40]. According to [41,42], the PEEC method shows satisfactory accuracy in comparison with experimental results and with simulation results calculated by the MoM and the FDTD method considering several types of electrodes with different diameters and positions with respect to the soil surface and lengths even when more complex structures are analyzed. Several PEEC numerical approaches applied to transient phenomena can be found in literature [43][44][45][46][47]. However, the extension of application toward gas-insulated substation transient analysis remains a challenging task due to high complexity of metallic structures involved in the computational domain.
The aim of the following study was to provide a clear understanding of the transient ground potential rise, across the metallic structures located inside GIS building, during switching operations considering different gas-insulated substation configurations (one GIS configuration contains a certain number of GIS modules) in order to identify the suitable modeling technique achieving efficient computational efforts. Through the proposed analysis methods, a parametric analysis was employed in order to accurately quantify the transient response of the system when an additional GIS bus section is connected to the model, by adopting a PEEC approach. For numerical modeling and simulation, the XGSLab software (XGSLab 9.4.1.5 version, SINT Ingegneria, Bassano del Grappa, Italy) package was used [48].

•
Based on a rigorous state-of-the-art analysis, the authors wanted to highlight the drawbacks of the numerical methods applied toward GIS transient behavior analysis due to several important reasons. The paper aimed to address two major drawbacks of the electromagnetic modeling techniques: • The lack of a validated digital avatar of the three-dimensional GIS metallic ensemble and adjacent metallic structures; The computational effort limitations arising when the DE-based methods are employed by presenting the PEEC approach application on GIS transient behavior analysis thus allowing to consider of the entire substation during the simulation process.
Up to now, due to the fact that there is no available digital model representing the three-dimensional GIS enclosure, the transient response of the grounding system in the presence of the metallic enclosure, connected in several locations to the grid, hence creating closed current loops affecting the fault energy flowing throughout the substation, could not be assessed.

General Description of the System
The investigated system is a 110-kV substation with three-phase GIS modules placed in a dedicated building. The substation has a double bus configuration containing four GIS modules, further noted from BUS1 to BUS4, and a bus coupler (BC) arranged in a horizontal layout. The two-bus system of the double bus configuration is placed vertically, one on the top of the other (see Figure 1). and simulation, the XGSLab software (XGSLab 9.4.1.5 version, SINT Ingegneria, Bassano del Grappa, Italy) package was used [48]. Based on a rigorous state-of-the-art analysis, the authors wanted to highlight the drawbacks of the numerical methods applied toward GIS transient behavior analysis due to several important reasons. The paper aimed to address two major drawbacks of the electromagnetic modeling techniques: • The lack of a validated digital avatar of the three-dimensional GIS metallic ensemble and adjacent metallic structures; • The computational effort limitations arising when the DE-based methods are employed by presenting the PEEC approach application on GIS transient behavior analysis thus allowing to consider of the entire substation during the simulation process.
Up to now, due to the fact that there is no available digital model representing the threedimensional GIS enclosure, the transient response of the grounding system in the presence of the metallic enclosure, connected in several locations to the grid, hence creating closed current loops affecting the fault energy flowing throughout the substation, could not be assessed.

General Description of the System
The investigated system is a 110-kV substation with three-phase GIS modules placed in a dedicated building. The substation has a double bus configuration containing four GIS modules, further noted from BUS1 to BUS4, and a bus coupler (BC) arranged in a horizontal layout. The twobus system of the double bus configuration is placed vertically, one on the top of the other (see Figure  1). Inside the coaxial pipe, the phase conductors are located in a radial configuration, with respect to the geometric origin of the enclosure.
The grounding grid layout can be described as follows: A copper strip contour (40 × 5 mm) is located on the GIS platform surroundings to which each GIS module is connected through two copper grounding leads. On the inner wall of the GIS building an additional copper strip contour is installed at h = 0.3 m, which ensures the conduction paths toward the external grounding system. Inside the coaxial pipe, the phase conductors are located in a radial configuration, with respect to the geometric origin of the enclosure.
The grounding grid layout can be described as follows: A copper strip contour (40 × 5 mm) is located on the GIS platform surroundings to which each GIS module is connected through two copper grounding leads. On the inner wall of the GIS building an additional copper strip contour is installed at h = 0.3 m, which ensures the conduction paths toward the external grounding system. Three copper conductor contours are buried outside the GIS building at different depths with respect to the soil surface and different distances from the GIS building wall. Taking into account the very fast time distribution of the transient regime, the auxiliary equipment located near the GIS building (autotransformers, lightning protection systems, etc.) were not considered during the computational process. The geometrical characteristics and the implemented materials, adopted during the computational process, associated with the system components are presented in Table 1.

CAD Approach
The GIS metallic enclosure was modeled through six parallel aluminum bars interconnected in such a manner that the resultant geometry constitutes the GIS shell equivalent structure (hexagon-shape geometry) (see Figure 2). The thickness of the aluminum bars was considered to be equal with the distance between the inner and outer walls of the GIS enclosure, accounting also for the thickness of both metallic walls. The metallic flanges located at each GIS equipment enclosure ends were modeled through aluminum conductors with a similar diameter as the aluminum bars, in order to ensure the galvanic connection between every two consecutive bars. Three copper conductor contours are buried outside the GIS building at different depths with respect to the soil surface and different distances from the GIS building wall. Taking into account the very fast time distribution of the transient regime, the auxiliary equipment located near the GIS building (autotransformers, lightning protection systems, etc.) were not considered during the computational process. The geometrical characteristics and the implemented materials, adopted during the computational process, associated with the system components are presented in Table 1.

CAD Approach
The GIS metallic enclosure was modeled through six parallel aluminum bars interconnected in such a manner that the resultant geometry constitutes the GIS shell equivalent structure (hexagonshape geometry) (see Figure 2). The thickness of the aluminum bars was considered to be equal with the distance between the inner and outer walls of the GIS enclosure, accounting also for the thickness of both metallic walls. The metallic flanges located at each GIS equipment enclosure ends were modeled through aluminum conductors with a similar diameter as the aluminum bars, in order to ensure the galvanic connection between every two consecutive bars. To understand and to endorse the adopted hexagon-based geometry approach, three simplified CAD geometries (pentagon, hexagon, and octagon layouts) were designed and tested in similar modeling conditions for a single GIS bus section. The endorsement and validation of the geometric approach were performed considering lightning surge scenario, considering a voltage breakdown fault. The purpose of the study was to quantify how the different number of parallel aluminum bars representing the solid metallic pipe will impact the transient response of the enclosure during voltage breakdown fault, taking into account very-high-frequency transients.
From a geometric point of view, each CAD model was built following similar design procedures. The radius of geometric circle of the enclosure (located at a half distance between the inner and the To understand and to endorse the adopted hexagon-based geometry approach, three simplified CAD geometries (pentagon, hexagon, and octagon layouts) were designed and tested in similar modeling conditions for a single GIS bus section. The endorsement and validation of the geometric approach were performed considering lightning surge scenario, considering a voltage breakdown fault. The purpose of the study was to quantify how the different number of parallel aluminum bars representing the solid metallic pipe will impact the transient response of the enclosure during voltage breakdown fault, taking into account very-high-frequency transients.
From a geometric point of view, each CAD model was built following similar design procedures. The radius of geometric circle of the enclosure (located at a half distance between the inner and the Energies 2020, 13, 6231 6 of 19 outer GIS walls), r = 0.248 m, was taken as reference, with all three investigated geometries being inscribed in a circle with the corresponding circle.
The average errors computed analyzing the numerical values of maximum and minimum amplitudes measured for the first through fourth periods of the time domain waveform were between acceptable limits, between 2-5%. Due to the fact that the spatial distribution of the aluminum bars surrounding the phase conductor was slightly different for each particular model, dissimilarities between numerical values were to be expected. However, the neglectable differences between recorded TGPR values considering different models were obtained as a result of conductive coupling, which dominates the transient response of the system during voltage breakdown fault. According to graphical as well as numerical results (see Figure 3), there were negligible differences between the primary parameters of the TGPR when different enclosures were considered during the computational process. Due to geometric symmetry reasons as well as from a computational time point of view, the hexagon-based geometry was developed and further used for the real GIS model analysis.
Energies 2020, 13, x FOR PEER REVIEW 6 of 19 outer GIS walls), r = 0.248 m, was taken as reference, with all three investigated geometries being inscribed in a circle with the corresponding circle. The average errors computed analyzing the numerical values of maximum and minimum amplitudes measured for the first through fourth periods of the time domain waveform were between acceptable limits, between 2-5%. Due to the fact that the spatial distribution of the aluminum bars surrounding the phase conductor was slightly different for each particular model, dissimilarities between numerical values were to be expected. However, the neglectable differences between recorded TGPR values considering different models were obtained as a result of conductive coupling, which dominates the transient response of the system during voltage breakdown fault. According to graphical as well as numerical results (see Figure 3), there were negligible differences between the primary parameters of the TGPR when different enclosures were considered during the computational process. Due to geometric symmetry reasons as well as from a computational time point of view, the hexagon-based geometry was developed and further used for the real GIS model analysis.     outer GIS walls), r = 0.248 m, was taken as reference, with all three investigated geometries being inscribed in a circle with the corresponding circle. The average errors computed analyzing the numerical values of maximum and minimum amplitudes measured for the first through fourth periods of the time domain waveform were between acceptable limits, between 2-5%. Due to the fact that the spatial distribution of the aluminum bars surrounding the phase conductor was slightly different for each particular model, dissimilarities between numerical values were to be expected. However, the neglectable differences between recorded TGPR values considering different models were obtained as a result of conductive coupling, which dominates the transient response of the system during voltage breakdown fault. According to graphical as well as numerical results (see Figure 3), there were negligible differences between the primary parameters of the TGPR when different enclosures were considered during the computational process. Due to geometric symmetry reasons as well as from a computational time point of view, the hexagon-based geometry was developed and further used for the real GIS model analysis.

Adopted Electromagnetic Field Theory Approach: Partial Equivalent Element Circuit Method
The partial element equivalent circuit (PEEC) method is derived from Maxwell's equations and provides a full-wave solution to them [41]. The PEEC method interprets Maxwell's equations to a circuit domain. The theoretical derivation of the PEEC method for the thin wire structures proposed by Yutthagowith and Ametani in [42] starts from a total electric field on a wire surface: Whole system equations in the frequency domain can be written in a matrix form corresponding to a modified nodal analysis (MNA) formulation, as shown in Equation (2): where A is an incident matrix that expresses the cell connectivity, R is a matrix of series resistances of current cells, L is a matrix of partial inductances of current cells including the retardation effect, P is a matrix of partial potential coefficients of the potential cells including the retardation effect, Y is a vector of potentials on the potential cells, I is a vector of currents along with current cells, U S is a vector of voltage sources, I S is a vector of external current sources, and Y a is an additional admittance matrix of linear and nonlinear elements [25]. The purpose of the computer-based solvers is to accurately expend analytical methods elaborated for basic structures to real-case scenarios considering metallic towers, different geometries of the grounding grid, gas-insulated substation metallic enclosure, and the electromagnetic interaction between such structures. Although, initially, the PEEC numerical approach was intended to be applied on inductivity analysis across integrated circuit boards, improvements and modifications have been made in recent years in order to be successfully applied on transient regime analysis. As was presented in previous sections, PEEC approach has several advantages compared to DE form-based methods, which require the discretization of the entire computational domain's volume, thus imposing limitations regarding the analysis of complex geometries. Moreover, the outcome provided by PEEC approach does not require post-processing procedures in order to obtain the potential and currents along a circuit as FDTD, MoM, and FEM methods. Regarding the application toward gas-insulated substation transient behavior analysis, there are barely a few papers throughout public literature focused on this topic.
The substation is energized by two 220/110-kV autotransformers, which are modeled through equivalent impedances, Z = 0.1 Ω. Pure resistive loads were assumed during the computational process, modeled as transversal impedances to Earth connected at each particular phase conductor Z = 1 Ω. The VFTO source was modeled through an ideal electromotive force (EMF) generator situated at the tripped disconnector switch location. Beside the transient source, each metallic element became an individual electromagnetic source, which generates a particular electric and magnetic field, affecting the metallic structures located in nearby surroundings, the electromagnetic chain reaction. The transient waveform is described by a double exponential function, presented in [12] with the following parameters α = 2.31049 · 10 5 s −1 , β = 8.17350372 · 10 9 s −1 , and V m = 1 p.u. in order to obtain per unit results considering fault severity. The following relation describes the time domain distribution of the voltage source in the XGSLab software package: where τ 1 = 1 β and τ 2 = 1 α with: k = correction factor, V m = maximum amplitude of the voltage waveform, τ 1 = time to peak parameter, and τ 2 = time to half value.
The proposed analysis was performed assuming a phase-to-enclosure fault inside of BUS4 (see Figure 1 for BUS4 location inside GIS).

Assessment Methodology
To correctly quantify the contribution of each particular metallic section at the holistic transient behavior of the substation, a parametric analysis was employed as follows: The first stage of the computational process will contain the initial GIS configuration (5 GIS bus sections, see Figure 1), taking into account phase-to-enclosure fault inside BUS4 enclosure, modeled through a galvanic connection between the phase conductor and the inner wall of the enclosure. The results obtained in this first analysis stage will represent the baseline for the investigation carried out in this paper.
With each future analysis step an additional bus section will be removed from the model until the single GIS bus section configuration is achieved (see  The proposed analysis was performed assuming a phase-to-enclosure fault inside of BUS4 (see Figure 1 for BUS4 location inside GIS).

Assessment Methodology
To correctly quantify the contribution of each particular metallic section at the holistic transient behavior of the substation, a parametric analysis was employed as follows: The first stage of the computational process will contain the initial GIS configuration (5 GIS bus sections, see Figure 1), taking into account phase-to-enclosure fault inside BUS4 enclosure, modeled through a galvanic connection between the phase conductor and the inner wall of the enclosure. The results obtained in this first analysis stage will represent the baseline for the investigation carried out in this paper.
With each future analysis step an additional bus section will be removed from the model until the single GIS bus section configuration is achieved (see Figures 5 and 6 for analysis of step 2 and step 4). Minimizing the number of elements contained by the model will reduce the computational effort bringing improvements of the required computational time. With each GIS bus section extracted from the computational domain, the number of elements contained by the model decrease by approximately 16% (considering that the complete configuration contained five sections and the grounding grid structure remained constant during the simulations).  Through overlapped graphical representations and numerical comparisons of the TGPR waveform computed for different GIS configurations, the transient response of the substation during voltage breakdown fault is assessed. This analysis will determine if the computed values of transient ground potential rise generated throughout the substation during phase-to-enclosure fault considering a single GIS bus section (the last stage of the investigation, see Figure 7) provide the worst-case scenario, within acceptable deviation range from the full GIS configuration model. As a result, a suitable modeling technique will be achieved and presented, and, therefore, the computational domain could be reduced for further analysis and investigations. Moreover, the transient response of the grounding grid in the presence of the metallic enclosure will be assessed and discussed. Due to the fact that there was no available 3D model representing the entire GIS metallic ensemble, quantifying the impact of the Energies 2020, 13, 6231 9 of 19 enclosure on the grounding grid performance when very-high-frequency transients flowed throughout the system was an impossible task. The response of the grid was expected not to be uniform across the substation, although it is based on the approximative symmetrical geometric shape of the earthing system due to the presence of the metallic enclosure. It is well known that during the very fast transient regime the effective area of the grounding grid is significantly reduced and, therefore, the grounding grid subsystems responsible for the fault energy clearance will be identified.    Through overlapped graphical representations and numerical comparisons of the TGPR waveform computed for different GIS configurations, the transient response of the substation during voltage breakdown fault is assessed. This analysis will determine if the computed values of transient ground potential rise generated throughout the substation during phase-to-enclosure fault considering a single GIS bus section (the last stage of the investigation, see Figure 7) provide the worst-case scenario, within acceptable deviation range from the full GIS configuration model. As a result, a suitable modeling technique will be achieved and presented, and, therefore, the computational domain could be reduced for further analysis and investigations. Moreover, the transient response of the grounding grid in the presence of the metallic enclosure will be assessed and discussed. Due to the fact that there was no available 3D model representing the entire GIS metallic ensemble, quantifying the impact of the enclosure on the grounding grid performance when very-high-frequency transients flowed throughout the system was an impossible task. The response of the grid was expected not to be uniform across the substation, although it is based on the approximative symmetrical geometric shape of the earthing system due to the presence of the metallic enclosure. It is well known that during the very fast transient regime the effective area of the grounding grid is significantly reduced and, therefore, the grounding grid subsystems responsible for the fault energy clearance will be identified.   Figure 7, represent the established analysis locations across the substations, as follows: upper side of the BUS4 (faulted bus) enclosure, grounding lead that connects BUS4 enclosure with the earthing system, copper strip located on the GIS platform, copper strip located on the inner wall of the GIS building, and vertical rod.
Quantifying the total discharge current by the enclosure toward the grounding grid considering such complex electromagnetic interactions within several metallic structures requires a detailed analysis of the contribution of each particular bus section. The entire structure of the grounding grid was considered during the computational process. Hence, the impact of adding or subtracting particular GIS bus sections from the model was analyzed.   Figure 7, represent the established analysis locations across the substations, as follows: upper side of the BUS4 (faulted bus) enclosure, grounding lead that connects BUS4 enclosure with the earthing system, copper strip located on the GIS platform, copper strip located on the inner wall of the GIS building, and vertical rod.
Quantifying the total discharge current by the enclosure toward the grounding grid considering such complex electromagnetic interactions within several metallic structures requires a detailed analysis of the contribution of each particular bus section. The entire structure of the grounding grid was considered during the computational process. Hence, the impact of adding or subtracting particular GIS bus sections from the model was analyzed.

Simulations and Results
The simulation was performed using an I7-7700 CPU, 3.60 GHz, personal computer with 16.0 GB of RAM and required 4 h of computational time for an 800-µs simulation period for the analyzed five GIS buses' configuration.
In the following section, several simplifications of the computational domain are proposed, tested, and performed in order to try to reduce computation time. Moreover, reducing the number of elements contained by the model will simplify the mathematical formulation associated with the electromagnetic field problem under study. When GIS arrangements contain four, three, two, and single modules then the computational effort has been reduced accordingly.
It has to be mentioned that when large matrix systems describing the model are conceived by the electromagnetic algorithm implemented in XGSLab software package, a well-known problem, known as an ill-conditioned system, might arise due to the low-frequency breakdown numerical instability. An ill-conditioned problem defines a condition of a small change in the inputs (the right-hand side of the equations' system) that leads to a large change in the output without real correlation with the physical phenomenon. However, using an appropriate numerical solver and suitable mathematical techniques, an accurate solution can be provided even for matrix systems containing a large number of elements and unknowns [49]. When low frequencies were considered during the simulation process, an artificial phase shift of the transient electromagnetic wave was observed, caused by the numerical instabilities associated with the low-frequency breakdown [50]. In order to overcome the low-frequency breakdown challenge, frequencies between 0-1000 Hz were not considered during the computational process. Considering the very fast nature of the physical phenomenon, the output of the method was not affected by the low-frequency spectrum elimination. Figure 8 shows characteristics of TGPR on the grounding lead associated with BUS4 under influence of faulted bus enclosure with the earthing system considering five, four, and one GIS bus section's(s') configuration. The analysis was performed on the grounding lead that connects the faulted bus enclosure with the earthing system for different gas-insulated substation configurations.

Simulations and Results
The simulation was performed using an I7-7700 CPU, 3.60 GHz, personal computer with 16.0 GB of RAM and required 4 h of computational time for an 800-μs simulation period for the analyzed five GIS buses' configuration.
In the following section, several simplifications of the computational domain are proposed, tested, and performed in order to try to reduce computation time. Moreover, reducing the number of elements contained by the model will simplify the mathematical formulation associated with the electromagnetic field problem under study. When GIS arrangements contain four, three, two, and single modules then the computational effort has been reduced accordingly.
It has to be mentioned that when large matrix systems describing the model are conceived by the electromagnetic algorithm implemented in XGSLab software package, a well-known problem, known as an ill-conditioned system, might arise due to the low-frequency breakdown numerical instability. An ill-conditioned problem defines a condition of a small change in the inputs (the righthand side of the equations' system) that leads to a large change in the output without real correlation with the physical phenomenon. However, using an appropriate numerical solver and suitable mathematical techniques, an accurate solution can be provided even for matrix systems containing a large number of elements and unknowns [49]. When low frequencies were considered during the simulation process, an artificial phase shift of the transient electromagnetic wave was observed, caused by the numerical instabilities associated with the low-frequency breakdown [50]. In order to overcome the low-frequency breakdown challenge, frequencies between 0-1000 Hz were not considered during the computational process. Considering the very fast nature of the physical phenomenon, the output of the method was not affected by the low-frequency spectrum elimination. Figure 8 shows characteristics of TGPR on the grounding lead associated with BUS4 under influence of faulted bus enclosure with the earthing system considering five, four, and one GIS bus section's(s') configuration. The analysis was performed on the grounding lead that connects the faulted bus enclosure with the earthing system for different gas-insulated substation configurations.  Considering additional GIS bus sections into the computational domain is equivalent to generating multiple parallelism conditions between each pair of metallic elements contained by the model (aluminum bar-aluminum bar, aluminum bar-copper strip) through which electromagnetic couplings are developed. While single GIS bus module configuration (only the faulted bus section) was analyzed, the transient response computed at the level of grounding lead was more severe (approximately double), due to the fact that the fault energy was not cleared throughout a complex metallic structure as in the case of four and five GIS bus sections' configuration. It can be noted that the computed voltage waveforms were characterized by similar harmonics when four and five GIS modules were considered during the computational domain.
When the analysis was performed on the upper side of BUS4 metallic enclosure, similar behavior to that of the TGPR was observed, as in the previous case: Gradually extracting additional GIS bus section from the model caused an amplifying effect on the oscillatory character of the voltage waveform combined with the reduction of the maximum amplitude (see also Figures 9 and 10). A steeper character 0.045 (time to peak) µs can be observed when a single bus is considered in comparison with 0.14 µs when four and five buses are considered, respectively. Similar harmonics and time-to-peak parameters describing the transient voltage waveforms computed during the computational process considering five and four GIS bus sections, respectively, could be observed (see Figure 7).
Energies 2020, 13, x FOR PEER REVIEW 11 of 19 Considering additional GIS bus sections into the computational domain is equivalent to generating multiple parallelism conditions between each pair of metallic elements contained by the model (aluminum bar-aluminum bar, aluminum bar-copper strip) through which electromagnetic couplings are developed. While single GIS bus module configuration (only the faulted bus section) was analyzed, the transient response computed at the level of grounding lead was more severe (approximately double), due to the fact that the fault energy was not cleared throughout a complex metallic structure as in the case of four and five GIS bus sections' configuration. It can be noted that the computed voltage waveforms were characterized by similar harmonics when four and five GIS modules were considered during the computational domain.
When the analysis was performed on the upper side of BUS4 metallic enclosure, similar behavior to that of the TGPR was observed, as in the previous case: Gradually extracting additional GIS bus section from the model caused an amplifying effect on the oscillatory character of the voltage waveform combined with the reduction of the maximum amplitude (see also Figures 9 and 10). A steeper character 0.045 (time to peak) μs can be observed when a single bus is considered in comparison with 0.14 μs when four and five buses are considered, respectively. Similar harmonics and time-to-peak parameters describing the transient voltage waveforms computed during the computational process considering five and four GIS bus sections, respectively, could be observed (see Figure 7).   Besides the parallelism conditions generated by multiple GIS sections contained in the computational domain, the proportions of equivalent inductance, capacitance, and resistance into the circuit were modified. From the harmonic point of view, the multiple parallel GIS modules behaved as a filter (see Figures 8,9 and 11).   The amplification effect influencing the oscillatory character of the transient overvoltage waveform provided by gradually extracting bus sections from the model could be observed also across the copper strip ring located on GIS concrete platform surroundings ( Figure 11). However, as the distance from the transient source increased, as well as from the aluminum bars' configuration,    The amplification effect influencing the oscillatory character of the transient overvoltage waveform provided by gradually extracting bus sections from the model could be observed also across the copper strip ring located on GIS concrete platform surroundings ( Figure 11). However, as the distance from the transient source increased, as well as from the aluminum bars' configuration, Figure 11. Transient ground potential rise across copper strip located on GIS platform, considering different GIS buses' configuration.
The amplification effect influencing the oscillatory character of the transient overvoltage waveform provided by gradually extracting bus sections from the model could be observed also across the copper strip ring located on GIS concrete platform surroundings ( Figure 11). However, as the distance from the transient source increased, as well as from the aluminum bars' configuration, the impact of electromagnetic couplings developed between particular metallic enclosure elements did not affect in a similar manner the maximum TGPR values recorded as in previous locations.
Let us assume that the first significant frequency period can be determined within the first 0.2 µs (nonperiodic signal) in all the cases presented in Figure 11. The time-to-peak parameters related to TGPR waveforms computed while five, four, and one GIS bus section(s) were 0.0525 µs for first two cases and 0.0775 µs, respectively, however, with different polarities. The multiple reflections of the electromagnetic wave travelling within aluminum bars' configuration, considering multiple GIS bus sections contained by the model, generated relatively steeper voltage waveform in comparison with a single GIS bus configuration. Figure 12 illustrates the transient overvoltage computed on the inner wall of the GIS building, analysis location, highlighted in Figure 5, considering different gas-insulated substation configurations.
Energies 2020, 13, x FOR PEER REVIEW 14 of 20 the impact of electromagnetic couplings developed between particular metallic enclosure elements did not affect in a similar manner the maximum TGPR values recorded as in previous locations. Let us assume that the first significant frequency period can be determined within the first 0.2 μs (nonperiodic signal) in all the cases presented in Figure 11. The time-to-peak parameters related to TGPR waveforms computed while five, four, and one GIS bus section(s) were 0.0525 μs for first two cases and 0.0775 μs, respectively, however, with different polarities. The multiple reflections of the electromagnetic wave travelling within aluminum bars' configuration, considering multiple GIS bus sections contained by the model, generated relatively steeper voltage waveform in comparison with a single GIS bus configuration. Figure 12 illustrates the transient overvoltage computed on the inner wall of the GIS building, analysis location, highlighted in Figure 5, considering different gas-insulated substation configurations. Considering the five GIS bus sections' configuration, the maximum amplitude describing TGPR computed across copper strip located on the inner wall of GIS building was attenuated with 87%, if compared with the maximum amplitude computed throughout the substation, which was considered as a reference value (on the metallic enclosure of BUS4).
When the transient electromagnetic wave reached the grounding grid, conductors located outside GIS building, vertical rods, and copper conductor, the oscillatory character was fully damped, as depicted in Figure 13, regardless of the GIS configuration considered during the computational process. Similar maximum amplitudes of the transient ground potential rise were computed on the vertical rod in comparison with those obtained at the previous analysis location; hence, the efficiency of the grounding grid decreased when the distance from the transient source increased. Comparing the TGPR computed, taking into account different GIS configurations, the worst-case scenario was obtained when five GIS buses were included in the computational domain and similar transient response of the vertical rods was computed when single and four GIS buses' configurations were considered. However, in order to achieve an accurate assessment regarding the holistic behavior of the substation, the number of GIS sections in the model must be limited in the computational process. Considering the five GIS bus sections' configuration, the maximum amplitude describing TGPR computed across copper strip located on the inner wall of GIS building was attenuated with 87%, if compared with the maximum amplitude computed throughout the substation, which was considered as a reference value (on the metallic enclosure of BUS4).
When the transient electromagnetic wave reached the grounding grid, conductors located outside GIS building, vertical rods, and copper conductor, the oscillatory character was fully damped, as depicted in Figure 13, regardless of the GIS configuration considered during the computational process. Similar maximum amplitudes of the transient ground potential rise were computed on the vertical rod in comparison with those obtained at the previous analysis location; hence, the efficiency of the grounding grid decreased when the distance from the transient source increased. Comparing the TGPR computed, taking into account different GIS configurations, the worst-case scenario was obtained when five GIS buses were included in the computational domain and similar transient response of the vertical rods was computed when single and four GIS buses' configurations were considered. However, in order to achieve an accurate assessment regarding the holistic behavior of the substation, the number of GIS sections in the model must be limited in the computational process. Energies 2020, 13, x FOR PEER REVIEW 14 of 19 Figure 13. TGPR, one particular vertical rod, comparison considering several GIS configurations.

Discussions
In general, if the assessment of safety procedures is successfully achieved for the worst-case scenarios, as a conservative measure, it can be said that the substation is fully protected. However, this could result in an overdesign of the required protection system and additional unnecessary financial costs. In order to achieve a safety level of the substation, the main characteristics of the protection system must overcome the transient overvoltage generated across the substation that is described by frequency, amplitude, and time-to-peak parameter. Designing a suitable safety protection system requires that the computational domain must contain the on-site physical configuration.
It is important to mention that the TGPR waveform illustrated throughout the current paper contained incident, reflected, and refracted components. The XGSLab software package computed the electromagnetic interaction between each pair of two metallic elements contained by the model. Therefore, the outcome of the method contained the combined effects of inductive, capacitive, and conductive couplings, fully taken into account.
Based on the conclusions from the state-of-the-art section (Introduction), a 3D model representing the entire metallic ensemble was required to accurately assess the holistic transient response of a gas-insulated substation and adjacent metallic structures during voltage breakdown fault, considering the very fast transient regime. The main concept of the proposed modeling technique was to represent the metallic enclosure (coaxial pipe) through several parallel aluminum bars. A first step of the proposed analysis method was to quantify how different numbers of parallel aluminum bars representing the solid metallic pipe will impact the transient response of the enclosure during voltage breakdown fault, taking into account very-high-frequency transients. According to the computed results (see Figure 3) there were negligible differences between the primary parameters of the TGPR when different enclosures were considered during the computational process, considering certain operating conditions.
Keeping in mind the hexagonal geometric model proposed, identified as most suitable for large GIS modeling, the transient ground potential rise was computed at several locations across the substations. According to the graphical and numerical results presented throughout the current paper, the spatial distribution of the TGPR across the grounding grid, especially inside GIS building, was not uniform, although the geometric arrangement of the earthing system was relatively symmetric with respect to the GIS platform, due to the presence of the metallic enclosure. As the distance from the transient source increased, as well as from the aluminum bars' configuration, the impact of electromagnetic couplings developed between particular metallic enclosures did not affect

Discussions
In general, if the assessment of safety procedures is successfully achieved for the worst-case scenarios, as a conservative measure, it can be said that the substation is fully protected. However, this could result in an overdesign of the required protection system and additional unnecessary financial costs. In order to achieve a safety level of the substation, the main characteristics of the protection system must overcome the transient overvoltage generated across the substation that is described by frequency, amplitude, and time-to-peak parameter. Designing a suitable safety protection system requires that the computational domain must contain the on-site physical configuration.
It is important to mention that the TGPR waveform illustrated throughout the current paper contained incident, reflected, and refracted components. The XGSLab software package computed the electromagnetic interaction between each pair of two metallic elements contained by the model. Therefore, the outcome of the method contained the combined effects of inductive, capacitive, and conductive couplings, fully taken into account.
Based on the conclusions from the state-of-the-art section (Introduction), a 3D model representing the entire metallic ensemble was required to accurately assess the holistic transient response of a gas-insulated substation and adjacent metallic structures during voltage breakdown fault, considering the very fast transient regime. The main concept of the proposed modeling technique was to represent the metallic enclosure (coaxial pipe) through several parallel aluminum bars. A first step of the proposed analysis method was to quantify how different numbers of parallel aluminum bars representing the solid metallic pipe will impact the transient response of the enclosure during voltage breakdown fault, taking into account very-high-frequency transients. According to the computed results (see Figure 3) there were negligible differences between the primary parameters of the TGPR when different enclosures were considered during the computational process, considering certain operating conditions.
Keeping in mind the hexagonal geometric model proposed, identified as most suitable for large GIS modeling, the transient ground potential rise was computed at several locations across the substations. According to the graphical and numerical results presented throughout the current paper, the spatial distribution of the TGPR across the grounding grid, especially inside GIS building, was not uniform, although the geometric arrangement of the earthing system was relatively symmetric with respect to the GIS platform, due to the presence of the metallic enclosure. As the distance from the transient source increased, as well as from the aluminum bars' configuration, the impact of electromagnetic couplings developed between particular metallic enclosures did not affect in a similar manner the maximum TGPR values recorded at various analysis locations. The maximum attenuation coefficients were recorded when the TGPR was analyzed starting from the upper side of the metallic enclosure up the GIS platform (69% when five-section configuration was considered, see Table 2). When the observation point moved toward the extremities of the grounding grid (namely, on the vertical rods) negligible TGPR attenuation rates were obtained. Therefore, it can be concluded that during very fast transients flowing throughout the substation the grounding grid effective area was located near the transient source, more precisely, on the GIS concrete platform. When the parametric study is employed several discussions can be initiated regarding the quantification of each GIS bus section into the overall transient response of the system. Through this type of analysis, the proper modeling technique is established. Moreover, through the number of GIS buses' variation during the simulations, the stability and modularity of the proposed modeling technique are achieved. The outcome of the method in all simulation conditions is strictly related to physical phenomenon rather than pure numerical dependencies. According to the TGPR comparison presented in a previous section (see Figures 9,11 and 12), the oscillatory character of the wave is damped within the traveling path due to resistive, inductive, and capacitive behavior of the metallic ensemble. Therefore, the frequency attenuation coefficient is directly proportional with the length of propagation path and equivalent with the number of elements contained by the model. According to the maximum overvoltage amplitudes presented in Table 2, when three, four, and five GIS sections are contained by the model, the area of the grounding grid clearing the transient energy is located inside the GIS building. Moreover, the copper strip situated beneath GIS metallic enclosure (copper strip ring and the grounding leads) will dominate the damping effect provided by the grounding grid during voltage breakdown (see Table 2).
When a single GIS bus configuration is considered during the computational process, the transient response of the system behaves differently, if compared with multiple sections' scenario. Similar TGPR maximum amplitudes are computed across the metallic enclosure and the grounding leads when a single bus section configuration is considered. Having multiple GIS sections in the model is equivalent to considering multiple existing electromagnetic couplings across the substation due to multiple parallelism conditions provided by the aluminum bars' ensemble, which greatly affect the overall transient response of the substation.
Different attenuation profiles of the transient electromagnetic wave phenomena are observed when different GIS configurations are considered during the computational process. According to the previous presented statements, it must be mentioned that, under no circumstance, should the simplification of the model in order to reduce the computational time be applied; whereas significant differences occurring between the topology of the computed time domain electromagnetic waveform in all considered configurations were observed (see Figures 9,11 and 12). Taking into account the different TGPR waveforms' topologies computed while single and multiple bus arrangements are considered, in order to design a suitable safety protection system, the computational domain must contain the on-site physical configuration.

Conclusions
Several important objectives were accomplished during the current study, which contributes to a better understanding of the holistic transient behavior of the gas-insulated substation.
The numerical instabilities due to the low-frequency breakdown were overpassed by neglecting the first frequencies' decades from the spectrum, between 0-1000 Hz during the computational process. By means of fast Fourier transform, the numerical algorithm embedded in the XGSLab software package computed the harmonic components of the transient source waveform, and the solution of the associated mathematical problem was obtained in the frequency domain. Applying the inverse fast Fourier transform, the time domain solution was computed. Considering the very fast nature of the physical phenomenon, the harmonic components with characteristic frequencies below 1000 Hz did not affect the overall transient behavior of the system under study.
A novel electromagnetic field-based model was developed, representing a real gas-insulated substation configuration. Based on the TGPR assessment at various locations across the substations, the proper modeling technique was established.
The effective area of the grounding grid during the voltage breakdown fault was identified. Based on the numerical and graphical results, a pattern regarding the fault energy flowing throughout the substation was established, which will contribute to a better understanding of the holistic transient behavior of a typical GIS substation.
The electromagnetic couplings developed inside the GIS building between the metallic enclosure and certain components of the grounding grid affected the behavior of the latter near the GIS platform.