Extended numerical simulation model for HVDC gas insulated systems considering gas humidity and gas convection

The growing global demand for electrical energy, combined with the construction of new power facilities in remote areas, has significantly increased the complexity of modern power systems. In addition, the necessity to meet environmental targets by using green energy, particularly hydropower plants, solar power systems and offshore wind farms, leads to a decentralization of the power grid and to larger transmission distances. Higher voltage (HV) levels are required to enable an efficient transmission of large electrical energy over long distances, in addition to direct current (DC) technology, since the reactive power losses under alternating current (AC) increase with transmission distance and beyond a certain line length AC becomes uneconomical (Zebouchi and Haddad, 2020), (Küchler, 2017). Therefore, High-voltage-direct-current (HVDC) gas insulated transmission lines (GIL) offer a reliable solution for these requirements.

HVDC GILs combine features such as high safety, reliability and a HV capacity compared to overhead lines or cables (Zebouchi and Haddad, 2020), (Zavattoni, 2014). Due to the excellent dielectric properties of the insulation gas, sulphur hexalfluoride (SF₆), which is compressed to pressures in the order of 1–6 bar to achieve a dielectric strength up to ten times higher than that of air (Kindersberger and Lederle, 2008), the GIL system requires considerably less space, making it suitable for offshore windplants or urban areas. In addition, GILs can be installed underground, resulting in a minimal visual impact on the environment. Therefore, the importance of HVDC GIL is expected to grow in the near future and will play a vital role in the transmission network of the power grid (CIGRE, 2021).

GIL consists of an inner conductor, an insulating gas, a grounded shell and a solid spacer, mainly made out of epoxy resin matrix filled with alumina (Al₂O₃) or silica (SiO₂), to sustain the system. The epoxy matrix is composed of a base resin with a hardener material and a catalyst. The most suitable variant is the aluminia-filled epoxy composite, due to its resistivity against SF₆ decomposition byproducts (Zebouchi and Haddad, 2020). There are different spacer geometries with different functions. Here, a HVDC GIL with a conical spacer is examined, which function is to stabilize the system and to separate the different gas chambers. The GIL geometry, including a conical spacer, is depicted in Figure 1. The interface between the gaseous and solid insulation represents the weak point of the system, since the different electrical properties of these materials may lead to space charge accumulation on the spacer surface, a phenomenon commonly observed under DC conditions. This slow charge accumulation may yield in locally high electric fields and partial discharges, endangering the safety of the system and generating electro-quasistatic conditions. Especially at the triple points between insulating gas, solid spacer and the electrodes (e.g. conductor or grounded enclosure), high electric fields occur. Hence, a precise understanding of the exact electric field distribution within the HVDC GIL is crucial for its optimal design. The charge accumulation phenomenon and its reduction by intelligent design of the HVDC GIL system is currently a subject of research (CIGRE, 2020), (CIGRE, 2021). Numerical simulations are a suitable method for investigating high voltage power equipment (Küchler, 2017). Therefore, an accurate simulation model of the HVDC GIL can be a powerful tool for design optimization. In most works, the electrical conductivity of the SF₆ gas is as assumed to be constant. However, measurements of the conduction phenomena in SF₆ gas reveal a complex conduction characteristic, with a nonlinear dependency on electric field, temperature, gas pressure and gas humidity (Hanna et al., 2016), (Zavattoni, 2014). These phenomena must be understood and accurately modeled to enable a reliable prediction of the space charge accumulation and high electric field regions in the system.

Thermal losses of the current in the conductor result in an increased temperature in the vicinity of the conductor, and therefore to a temperature gradient between conductor and grounded enclosure. Since the electrical conductivity of the SF₆ gas and the spacer material exhibits a nonlinear dependence on electric field and temperature, the electric field and temperature are intrinsically coupled through the conductivity of the insulating materials. Therefore, a coupled electro-thermal model incorporating a nonlinear electrical conductivity model for both the SF₆ gas and the spacer material must be developed. To be able to consider the influence of the temperature, the accurate temperature distribution has to be determined in the simulation model. Therefore convective heat transfer and heat radiation in the gaseous domain and heat conduction in the solid domain must be included in the model. This leads to a more complex model compared to most works, where heat conduction for the whole HVDC GIL system is considered or the temperature as a whole is neglected (Lavesson and Doiron, 2015), (Hayakawa et al., 2018a), due to the multiphysical characteristic of this approach (Vincenti and Kruger, 1975).

The aim of this study is to develop an accurate 3D simulation model for HVDC GIL that takes into account complex physical phenomena, based upon previous research (Hensel et al., 2022), (Hensel and Clemens, 2024). In (Hensel et al., 2022) an electrical conductivity model for SF₆ gas was developed for a 2D-axisymmetric simulation model with a nonlinear dependency on the electric field, the temperature and the gas pressure. This model was verified by ion-drift-diffusion simulations in Hensel and Clemens (2024), showing its accuracy with less computational cost. However, the simulation model in these works did not take into account the humidity content of the gas and only heat conduction phenomena was considered. Here, the electrical conductivity model is extended by also considering the gas humidity dependency and the simulation model is extended by developing a 3D model, which also takes heat convection and heat radiation into account, allowing to consider the effect of a vertical or horizontal arrangement of the GIL.

This work is organized as follows: after the introduction, the derivation of the electro-quasistatic equation is presented, followed by a description of the equations for calculating the temperature distribution and fluid dynamics within the system. In addition, since the temperature and the electric fields are coupled by the electrical conductivity, the nonlinear electrical conductivity model of the spacer material is introduced and the equation for the SF₆ gas’ electrical conductivity is derived, based on measurements. Subsequently, the results with the consideration of the gas humidity and its influence on the electrical conductivity of SF₆ -and consequently on the electric field distribution- are presented. Afterwards, the results including the effects of heat convection and thermal radiation on the temperature distribution and electric field are exhibited. Finally, the conclusions are provided.

The electro-quasistatic field formulation for the performed finite-element-method (FEM) simulations are subsequently derived (Jin, 2014), (Zienkiewicz et al., 2013). Electro-quasistatic fields are seen in applications where the electric field changes slowly and electromagnetic wave propagation and magnetic induction can be neglected. The system behaves nearly electrostatic with slow time variation. These conditions are given in HVDC gas insulated systems, since the wavelength is much larger than the dimension of the system and the electric field energy in high-voltage components is much larger than the magnetic field energy, therefore, inductive effects are negligible. When DC voltage is applied in the gas insulated system, charges accumulate in the insulating material, slowly changing the electric field over time. Therefore, under DC conditions solving the Laplace equation is not sufficient to determine the electric field distribution, apart from AC condition. In general, to calculate the electric field, the Poisson-equation, continuity equation and Ohm’s law are used:

where $J→$ [A/m²] is the current density, $ρ$ [As/m³] the space charge density, $ε0=8.854·10-12$As/(Vm) the vacuum permittivity, $εr$ the relative permittivity, $κ$ [S/m] the electric conductivity, T [K] the temperature, $φ$ [V] the electric potential and $E→$ [V/m] the electric field (De Lorenzi et al., 2009), (Luo et al., 2020). Under electro-quasistatic assumption of irrotational electric fields, $E→$ and $φ$ are in relation according to:

When (1), (3) and (4) are substituted into (2), the electro-quasistatic formulation is obtained:

Here, t [s] is the time (Balanis, 2012).

Since the electric conductivity of the insulating materials is decisive for the electric field and has a nonlinear dependency on temperature, the exact knowledge of the temperature distribution is crucial to determine the electric field distribution (Hayakawa et al., 2018a). In the solid domain (spacer material) the distribution of the temperature can be calculated by using the heat conduction equation:

where $κT,E→E→2$ describes the electric loss density in the insulation material, $λ$ [W/m$·$K] the thermal conductivity, $δ$ [kg/m³] the material density and $cp$ [J/kg$·$K] the specific heat capacity (Spitzner, 2013).

The determination of the temperature in the gaseous domain (SF₆) is more difficult, however. Unlike in the solid, where heat conduction can be considered, thermal convection is the dominant temperature transfer mechanism in the gas. Therefore, the Navier-Stokes equation needs to be solved to determine the flow velocity of the gas and hence the convection to identify the temperature distribution. The equation for thermal convection and the continuity equation are then solved, to determine the temperature in the gas. These equations are given by:

where u [m/s] is the flow velocity, $αp$ [1/K] the coefficient for the thermal expansion, $Tamb$ [K] ambient temperature, $P$ [MPa] the gas pressure, $μ$ [m²/s] coefficient of the kinetic viscosity, $I$ the identity matrix, $g$ [m/s²] the acceleration of gravity (Spitzner, 2013), (Vincenti and Kruger, 1975).

In addition, the effect of heat radiation between the electrodes is considered by solving the following equation:

where $σ$ [1/mol] is the Avogadro Boltzmann constant, $ο$ [W/(m²$·$K)] the radiation coefficient, n the refractive index, $Tamb$ [K] ambient temperature$.$

The electrical conductivity $κ$ is the decisive parameter for the electric field distribution under DC conditions. Therefore, the exact electrical conductivity value of the insulation materials must be known for the determination of an accurate electric field distribution. Measurements shows that for the epoxy spacer material, the electrical conductivity shows a nonlinear dependency on temperature and electric field and is described by:

where $kB=8.617·10-5$ eV/K is the Boltzmann constant, $WA=0.95$eV the activation energy and $κ0=8.5$ and $ϑ=0.08·10-6$ are constants (Hayakawa et al., 2018a).

The electrical conductivity of $SF6$ is commonly assumed to be constant in the majority of published models (De Lorenzi et al., 2009), (Hayakawa et al., 2018a), but is here also described by a nonlinear model. Measurements of the dark current behavior in SF₆ show a nonlinear dependency of the electrical conductivity on the electric field, temperature, gas pressure and the gas humidity (Zavattoni, 2014), (Hanna et al., 2016). Dark currents are small electrical currents that appear in insulating materials even without external ionizing mechanisms. They arise from natural ionization mechanisms such as background radiation. Although they are weak conduction currents, they may have an influence on the insulation and breakdown behavior. From these measurements a model with a nonlinear dependency on the electric field, temperature and gas pressure was derived in Hensel et al. (2022) and validated in Hensel and Clemens (2024). Here, the model is extended by also considering the nonlinear dependency on the humidity content of the gas. The current density in the SF₆ gas was measured during a variation of the electric field, the temperature, gas pressure and the gas humidity. These measurement values are depicted in Figures 2–5. The derivation of the model was conducted by empirically fitting a parametric function to the dark current measurements in an iterative process. The resulting current density function was divided by the electric field $E→$ to obtain the electrical conductivity according to (3). The derived parametric function, which describes the electrical conductivity of SF_6, is given by:

where: $κSF6=1·10-23, α = 40$⁠, $β = 4$⁠, $γ = -0.11$⁠, $Ex=21·106$⁠, $ς=12$⁠, $ϱ=0.8$⁠, $ϵ = 3.5$⁠, $Ey=16.5·106$⁠, $ι=13$⁠, $ζ = 2.1$⁠, $ν = 0.026$ and $ξ = 0.09$ are constants. The model allows to consider the effect different parameters, such as electric field and gas humidity, have on the electrical conductivity of the insulating gas and hence on the electric field distribution, which would be neglected if the electrical conductivity was set constant.

Resulting from various ionization processes such as background radiation and cosmic rays, a small dark current is flowing in the SF₆ gas. In addition, when a voltage is applied in the GIL, charge emission from the electrodes must be considered. Measurements from Zavattoni (2014) and Hanna et al. (2016) show that the electric conductivity of SF₆ has a nonlinear dependency on the electric field, temperature, gas humidity and gas pressure, which are depicted in Figures 2–5. From these measurements (12) is derived. The dependency of (12) on the parameters $E→$⁠, $P$⁠, $T$ and $RH$ is demonstrated in Figures 2–5, where the measured current density in SF₆ is shown and compared to (12). The electrical conductivity model is visualized by multiplying $E→$ to (12), to obtain the current density $J→$ according to (3).

The 3D geometry model of the examined HVDC GIL is shown in Figure 6. The simulations are performed with Comsol Multiphysics®, which uses the finite-element-method (Zienkiewicz et al., 2013). The height of the GIL is 1 m and the length is 2 m, as it can be seen in Figure 6. The radius of the spacer is 0.35 m and the conductor has a radius of 0.1 m. The HV conductor has a potential of $φ=320$kV and the shell is grounded (⁠$φ=0$kV). The relative permittivity of the epoxy resin and the insulating gas is assumed to be $εSpacer=4$ and $εSF6=1$⁠. Here, the permittivity values of the insulation materials are assumed to be constant, since under DC-conditions the electrical conductivity is the decisive factor for the electric field distribution. Since the frequency in HVDC GIL is effectively 0 Hz and fast transients do not occur, no dielectric dispersions appears and a static polarization dominates. Therefore, this assumption does not negatively affect the accuracy of the simulation model. In addition, under practically all operating conditions a highly constant relative permittivity can be assumed for SF₆ (Küchler, 2017). The temperature gradient between the conductor and the grounded shell is 40 K. At the truncated domain in z-direction the boundary condition is set as Neumann boundary condition, which means there are no normal electric field components to prevent artificial field distortions. In addition, the heat flux is zero in z-direction. The values of the boundaries represent typical conditions in a HVDC GIL (CIGRE, 2021), (Zebouchi and Haddad, 2020). The electrical conductivity models (11) and (12), for the spacer material and the SF₆gas, respectively, are applied.

The numerical simulation of the HDC GIL is a strongly coupled multiphysical problem that relies on bidirectional coupling between the electric, thermal and the flow fields. These quantities interact through Joule heating, buoyancy-driven convection and radiative heat transfer. The multiphysical problem exhibits strong nonlinear feedback mechanisms. The electric field determines Joule heating, which raises the temperature. The temperature increase influences the electrical conductivity of the insulation materials and the gas density, thereby influencing the electric field distribution and the buoyancy-driven flow again. The gas flow, in turn, redistributes the heat, while thermal radiation contributes additional surface heat exchange.

The coupled governing equations, discretized by using FEM, are solved by a segregated iterative approach, using the block Newton-method with a tolerance-based termination. Here, the problem is subdivided into several segregated steps, where each step represents the physics-based groups of dependent variables, which are solved sequentially. This strategy requires less memory and enables more numerical robustness, especially for 3D problems, compared to a fully coupled solver, where a single large system of equations is solved for all the unknowns within a single iteration. The termination criterion is reached when the relative tolerance exceeds the relative error. The model is discretized with an adaptive mesh. This solving approach represents a realistic prediction of the electro-thermal behavior of the system under HVDC operating conditions.

Here, the steady state condition is of interest, which is reached when the charge accumulation in the insulating material is terminated, resulting in a constant electric field distribution. The following examined results are at steady-state condition. For comparison, the simulation results of the electric field distribution in the HVDC GIL model assuming a constant electrical conductivity of $κSF6=1·10-18$ S/m for the SF₆ gas, as commonly applied in most studies Hayakawa, (2018a), are presented in Figure 7.

Figure 8 shows the electric field distribution in the HVDC GIL using the electrical conductivity model (12) for SF₆, assuming a gas pressure of $P=0.6 MPa$ and $RH=2$%, which is representative of dry gas conditions (Zavattoni, 2014) (left). It is notable that, in this case, only heat conduction is considered throughout the entire system, as commonly done in literature (Hayakawa et al., 2018a), (Hayakawa et al., 2018b). Simulation results with consideration of heat convection and radiation are presented subsequently, for a better comparison.

The electric field shows an inversion in the vicinity of the spacer, due to the charge accumulation under DC condition and the temperature gradient between the grounded shell and the HV conductor. This is due to the temperature dependency of the electrical conductivity of the insulation materials. Since the temperature at the conductor is increased, the electrical conductivity of the spacer material is also increased in the vicinity of the conductor, which leads to a decrease of the electric field in this environment. Therefore, higher electric fields occur at the grounded enclosure, especially at the triple point. It is notable that the electric field is concentrated at the upper side of the spacer. Compared to Figure 7, where the electrical conductivity of SF₆ is assumed to be constant, influences of parameters such as temperature, gas pressure and gas humidity are neglected and consequently, the simulation results differ from those shown in Figure 8. There, higher electric fields occur, with a maximum value of 4.03 kV/mm, compared to the 3.88 kV/mm in Figure 8 (left). However, the largest difference can be seen in the electric field distribution, where considerably higher electric fields occur in the vicinity of the conductor, compared to Figure 8. This behavior is attributed to the higher electric field strength near the conductor, which, by applying (12), results in an increased electrical conductivity in this region and in reduced local electric field strengths. This demonstrates that neglecting the nonlinear behavior of the electrical conductivity of SF₆ leads to less accurate simulation results.

This effect becomes even more distinct when parameters such as the gas humidity exceed typical levels. Measurements on humid SF₆ gas reported by Zavattoni (2014) indicate that a representative humidity value is approximately 30%, which is applied on the right-hand side in Figure 8. In comparison with dry gas conditions (3.88 kV/mm), the maximum electric field differs by approximately 35% and is with 2.51 kV/mm decreased. This behavior is attributed to the increased electrical conductivity of SF₆, resulting from its higher humidity content.

In the following, simulation results that additionally account for heat convection in the gas and thermal radiation are discussed-by solving equations (7)–(10) are discussed. Even though heat convection is the dominant temperature transfer mechanism in the gas compartment of HVDC GILs, heat radiation is still considered, for the reason that high accuracy is desired for the simulation model. Since convection and the gas flow dynamics are influenced by gravity, it is necessary to distinguish between a horizontally oriented GIL and a vertically arranged GIL.

Figure 9 depicts the gas flow velocity for a horizontally arranged GIL, where gravity is applied in negative y-direction, and for a vertical arranged GIL, where gravity is applied in negative z-direction. Even though GIL chambers are mechanically separated by cone-type spacers, a gas velocity can still occur. The spacer provides electrical insulation, but due to several physical mechanisms, a fluid flow still appears. Under high electric fields, ions drift in the field and due to momentum transfer on gas molecules, a gas flow results. Also due to the temperature gradient between grounded enclosure and the conductor a convective circulation occurs, accompanied by pressure fluctuations and surface guided flows along the spacer, resulting in gas motion. In Figure 9 it can be seen that there is a distinct difference in the flow velocity of the gas, depending on a horizontally or vertically arranged GIL. In case of the horizontally positioning there is a higher flow velocity seen in the upper part between spacer and conductor, resulting from a different direction of the gravity force, compared to the vertical GIL. This phenomenon commonly occurs in practical applications, where the GIL may be installed horizontally on the ground, buried underground or oriented vertically, depending on the environmental conditions. Therefore, it is important to account for fluid dynamics and the associated heat convection to capture these effects (Hu et al., 2023), (Neidhart et al., 2020).

The corresponding temperature distributions in both cases can be observed in Figure 10. According to Figure 10, for a horizontal GIL, the temperature is increased in the upper parts between spacer and conductor. In contrast, for a vertical GIL, the temperature distribution is more homogenous, with higher values concentrated near the conductor.

Due to the nonlinear dependence of the insulation materials’ electrical conductivities on temperature, and because the electrical conductivity directly affects the electric field, the electric field distribution should be investigated for the different scenarios. Figure 11 shows the electric field distribution for dry gas (RH = 2%) including heat convection and thermal radiation, for a horizontally oriented GIL (left) and a vertically arranged GIL (right).

Compared to Figure 8, where heat convection and thermal radiation are not considered, the electric field exhibits a different distribution and a different maximum value. The electric field is still inversed, however, the electric field is more increased in the vicinity of the conductor and the maximum value with 2.89 kV/mm (horizontal) and 2.34 kV/mm (vertical) is considerably lower compared to the 3.88 kV/mm from Figure 8. This further demonstrates the importance of accurately representing the physical heat transfer mechanisms in HVDC GIL simulation models. Furthermore, the maximum electric field also exhibits a clear dependence on the GIL orientation.

The electrical conductivity model (12) for SF₆ enables the investigation of the influence of gas humidity on the electric field distribution. Accordingly, the electric field distribution in humid SF₆ (RH = 30%) is analyzed using the HVDC GIL simulation model, which accounts for gas convection and thermal radiation, and is compared to the results shown in Figure 8. Figure 12 presents the simulation results for both scenarios. Similar to Figure 8, an increased humidity content in the SF₆ gas leads to a lower maximum electric field, due to the corresponding increase of the electrical conductivity of the SF₆ gas. The maximum electric field in the horizontally arranged GIL is decreased from 2.89 kV/mm to 2.13 kV/mm and in the vertical GIL from 2.34 to 1.83kV/mm. Similar to dry gas, the electric field distribution differs compared to Figure 8, where the temperature transfer mechanisms of heat convection and thermal radiation are neglected. Here, the maximum electric field reaches 2.51 kV/mm, which is considerably higher than the results shown in Figure 12. Interestingly, for a vertically arranged GIL, the electric field stress is lower compared to the horizontal configuration, for both dry and humid SF₆ gas. These observations further highlight the importance of incorporating such physical mechanisms in numerical simulations of HVDC GILs.

An important aspect of numerical computation is the solution time. Generally, a trade-off exists between the accuracy of a simulation model and the required solution time, where the goal is to achieve the desired precision without excessive computational cost. Table 1 presents the average computation times for the different scenarios examined in this study. The first model, which neglects fluid dynamics, heat convection and thermal radiation, requires an average of 3 min to obtain a solution. The model that incorporates these phenomena and solves the equations (7)–(10), exhibits varying solution times depending on the arrangement of the GIL. The vertical model requires 27 min and 30 s on average, whereas the horizontal model requires 37 min on average. Both models require significantly more (up to 12 times) computation times compared to the first scenario, however, they also allow the consideration of gas fluid dynamics and the orientation of the GIL.

In this work, an extended 3D numerical simulation model for a HVDC GIL was developed, incorporating complex physical phenomena. Unlike conventional models that assume a constant electrical conductivity value for the SF₆ gas, a nonlinear conductivity model, derived from dark current measurements in the gas, is applied. This model exhibits a nonlinear dependency on the electric field, temperature and gas pressure. Here, the model is extended by also considering a nonlinear dependency on the gas humidity content, as measurements indicate a direct effect on the conduction behavior of the gas. The simulation results demonstrate that variations in the humidity content have a significant influence on the electric field distribution and therefore should not be neglected in precise insulation design. The difference in the maximum electric field between dry and humid gas is up to 35%.

Furthermore, the model extends the thermal analysis by accounting not only for heat conduction as in most works, but also considers complex heat transfer mechanisms such as heat convection and radiation phenomena in the gas. This enables the coupling of fluid dynamic effects with temperature and electric field analysis, providing a multiphysical simulation model. The model also allows to consider the arrangement of the GIL, e.g. when the positioning is vertical or horizontal, leading to a different gas fluid dynamic. The simulation model results show different flow velocities, temperature behavior and electric field distributions depending on the configuration of the GIL, which is also seen in practical applications, confirming the importance of considering such effects in numerical simulation models.

The extension of the simulation model for HVDC GIL leads to a more complex model and therefore to larger solution times. The consideration of such effects leads to up to 12 times larger simulation times, compared to simulations where these effects are neglected. Overall, the proposed model represents a substantial improved simulation approach for HVDC GIL systems. It enables to consider effects such as gas humidity and the arrangement of the system and therefore a deeper investigation of design optimization.