Numerical Simulation of Temperature and Pressure Changes due to Partial Discharges in Spherical Cavities Within Solid Dielectrics at Different Ageing Conditions

Abstract

Partial Discharges (PD) behavior during ageing of the insulation systems exhibits variations that depend on changes in gas filling characteristics and surface condition. In this article, numerical simulations of temperature and pressure behavior in an air-filled spherical cavity within a homogenous solid dielectric material due to PD activity are presented. An Analytical-Finite Element Analysis simulation approach was implemented in MATLAB and results exhibit reasonable agreement with experimental measurements reported by other authors. Simulation results allow concluding that pressure changes are directly related to variations in the PD behavior. In addition, affectations to cavity surface due to temperature increments can be discarded.

1. Introduction

Partial Discharges (PD) are one of the main causes of dielectric breakdown in solid insulators and are also a symptom of ageing in high-voltage (HV) equipment [1,2]. Its simulation allows analyzing different conditions such as, waveform, ambient temperature and initial pressure, and geometrical configuration, among others. In addition, it allows to characterize the PD phenomena and determine the main parameters that affect the PD behavior under some specific conditions [3,4].

PD in cavities within solid dielectrics are modelled as streamers because they are self-sustained discharges and have greater magnitudes than other kinds of internal discharges [5]. There are different models for simulating the PD activity inside cavities within solid dielectrics and they differentiate in the way the electric field produced by the electric charge on the void surface, left by previous PD, is calculated [6]. PD models can be classified into [7]: analytical [8,9,10], three-capacitance [11] and finite-element analysis (FEA) [12,13] models.

In the three-capacitance models, the cavity is represented by an equivalent capacitance and PD are simulated reducing the magnitude of a resistance in parallel with that [14]. Analytical models use analytical expressions, that were theoretically determined using an electrostatic field approach [8], for calculating the electric field strength and the real and induced PD charge magnitude. On the other hand, FEA models use numerical methods for solving the partial differential equations (PDE) that describe the PD phenomena, and the real and induced charges are dynamically calculated integrating field expressions on HV electrodes boundaries [6].

Analytical models, in spite of their simplicity, are electrostatic field-based approaches and when homogeneous and symmetrical conditions are considered, such as cylinders or spheroids, the solutions obtained through analytical and FEA models are similar [15], and have the advantage that less computational time is required when simulations for a great number of periods of the AC-voltage signal are required such as in ageing analysis.

PD are inherently multi-physical phenomena where electrical, thermal and mechanical physics are involved in each PD event [16]. Processes such as electron generation rate and surface charge decay through surface conduction are temperature and pressure dependent [17,18,19]. On the other hand, temperature affects the ageing process in solid dielectrics mainly due to the acceleration of chemical ageing dynamics [1].

In this paper, the analytical and thermal models are presented in Section 2. In Section 3, PD in a spherical cavity are simulated using the analytical model at different ageing phases and results are compared with measured values reported by other authors. In addition, temperature and pressure variations due to PD events are calculated using a heat transfer FEA model, and the effect of temperature and pressure along ageing on the electron generation rate and in inception and extinction field magnitudes is analyzed and correlated with measured variables. A discussion on the simulation results is presented in Section 4. Finally, a conclusion is presented in Section 5.

2. PD Analytical and Temperature Models

In this study, an analytical-FEA approach is used for implementing simulations. With the analytical model, the electric phenomena related to PD events are simulated, and variables such as real and induced charge magnitudes are calculated [8]. On the other hand, temperature and pressure in gas and solid dielectrics due to PD activity are determined solving the heat transfer equation using a FEA model [20].

2.1. PD Stochastic Model

Through experimentation and analytical studies, it has been found that PD phenomena are stochastic processes where their properties are describable by time-dependent random variables [21]. The streamer development process depends on the achievement of the two following necessary conditions:

• The electric field strength magnitude inside the cavity is greater than the inception value;
• There is at least an electron for starting the avalanche.

The achievement of necessary conditions for PD occurrence outcome into the stochastic behavior of PD phenomena, determining characteristics such as inception delay, phase occurrence and the number of PD per cycle. The electric field strength inception magnitude is calculated using Equation (1) [9]:

$$E_{inc}={\left(E/p\right)}_{cr}p\left[1+\frac{B}{{\left(2pa\right)}^n}\right],$$

(1)

where $p$ is the pressure of gas within the cavity in Pa, $a$ is the cavity radius in m, ${\left(E/p\right)}_{cr}=24.2$ V·Pa⁻¹·m⁻¹, $B=8.6$ Pa^1/2·m^1/2 and $n=0.5$. The magnitude calculated using Equation (1) is the critical electric field strength magnitude for starting an avalanche.

Volume and surface emissions are the main mechanisms for the first electron generation rate [22]. Volume generation is related to radiative gas ionization and field detachment of electrons from negative ions and can be calculated as [23]:

$$N_{et}(t)=C_{rad}{\Phi }_{rad}{\left({\rho }_m/p\right)}_{0}p\left(\frac{4}{3}\pi a^3\right)\left(1-{\nu }^{-1/n}\right),$$

(2)

where $C_{rad}{\Phi }_{rad}=2\times {10}^6$ kg⁻¹·s⁻¹, is the reduced radiative cosmic and radioactive quantum flux density and ${\left({\rho }_m/p\right)}_0=1\times {10}^{-5}$ kg·m⁻³·Pa⁻¹, is the pressure reduced gas density.$\nu =U_{cav}(t)/U_{inc}$, $U_{cav}(t)$ is the voltage across the cavity center at instant $t$ and $U_{inc}$ is the inception voltage calculated with Equation (1).

The electron generation rate due to surface detrapping obeys the Richardson–Schottky [6] law and can be calculated using the following expression:

$$N_{dt}(t)=N_{dt0}\exp \left(-\frac{t-t_{PD}}{\tau }\right){\nu }_0\exp \left(-\frac{{\Phi }_{dt}-\sqrt{e{E}_{cav}(t)/(4\pi {\epsilon }_0)}}{k_{b}T}\right),$$

(3)

where $e$ is the elementary charge, ${\Phi }_{dt}$ is the effective detrapping work function, $k_b$ is the Boltzmann constant, $t-t_{PD}$ is the time elapsed since the latest PD, ${\nu }_0=1\times {10}^{14}$ Hz is the fundamental frequency of phonon, ${\epsilon }_0$ is the permittivity of vacuum and T the temperature in K. $N_{dt0}=\xi \left(q/e\right)$ and $\xi$ is a proportional factor, which describes the fraction of charge carriers that result in the creation of detrappable electrons. It depends on the relative polarity of electric field strength produced by charges on the surface related with the polarity of resultant electric field strength.

An electron is generated at the time interval $\left[t,t+\Delta t\right]$ with the probability $\left(N_{et}(t)+N_{dt}(t)\right)\Delta t$.

2.2. PD Analytical Model

In this model, it is considered that all the cavity volume is affected by a PD event and the electric field strength inside the cavity is uniform; for that reason, all calculations are referred to the cavity center. During a PD event, charges of different polarity are deployed on the void surface, which produces a dipole moment as shown in Figure 1.

In Figure 1, $s$ is the direction vector of dipole, pointing from negative to positive charges, $\sigma \left(\theta \right)$ is the surface charge density and $E_0$ is the electric field strength established by the HV electrodes inside the homogeneous dielectric material.

For ellipsoidal voids, the induced charge on the HV electrodes by the dipole in the cavity can be calculated as Equation (4) [9]:

$$q^{\prime }=-K\Omega \epsilon \left(E_{cavPD}-E_{ext}\right)\cdot \nabla {\lambda }_0,$$

(4)

where $K$ is a dimensionless factor dependent on the geometry of the cavity and for spherical cavities is $K=3$ [24], $\Omega$ is the cavity volume, ε is the permittivity of media, $\nabla {\lambda }_0$ is the electric field strength by unit of applied voltage in the homogeneous material without cavity and charges, and at the cavity center location. $E_{cavPD}$ is the electric field strength at the cavity center just before a PD event occurs. $E_{ext}$ is the electric field strength below there is no ionization and the streamer extinguishes, this can be calculated using Equation (5).

$$E_{ext}=\gamma {\left(E/p\right)}_{cr}p,$$

(5)

where $\gamma$ is a dimensionless factor depending on the polarity of the streamer. As average for positive and negative streamers, it is assumed as $\gamma =0.35$ [22].

The real charge on the void surface is proportional to the electric field change due to a PD event and can be calculated as Equation (6) [9,22]:

$$q={\epsilon }_0\pi a^2\left[1+{\epsilon }_r\left(K-1)\right)\right]\left(E_{cavPD}-E_{ext}\right),$$

(6)

where ${\epsilon }_r$ is the relative permittivity of the solid dielectric.

The charge magnitude on the cavity surface left by previous PD decay due to recombination, conduction and diffusion phenomena [25]. However, recombination due to surface currents is the dominant process [26]. The decay of charge as a function of time can be approximately described by the Ohm’s Law as in Equation (7) [23]:

$$-\frac{dq}{dt}=\left(\frac{\pi }{2}\right)k_s{E}_{cav}2a,$$

(7)

where $k_s$ is the cavity surface conductivity.

2.3. Temperature and Pressure Model for PD

During a PD event, the energy stored in the electric field changes abruptly, and using the Poynting’s theorem [27], it is concluded that the change in the accumulated energy is due to an energy conversion into heat which is transferred to cavity and solid material through heat conduction. The heat transfer can be simulated solving the following PDE in Equation (8) [20]:

$${\rho }_m{C}_p\left(\frac{\partial T}{\partial t}\right)-\nabla \cdot \left(k_T\nabla T\right)=Q,$$

(8)

where ${\rho }_m$ is the density of media, $C_p$ is the heat capacity of media and $k_T$ is the thermal conductivity of media. $Q$ is the heat source density and can be calculated as Equation (9):

$$Q=\frac{1}{2}\epsilon \left(\frac{{\left|E_{cavPD}\right|}^2-{\left|E_{ext}\right|}^2}{\Delta t_{PD}}\right),$$

(9)

where $\Delta t_{PD}$ is the time elapsed during the PD event. In this study, the PDE toolbox of MATLAB [28] is used for solving Equation (8), considering there is not a heat source density within the homogenous material and exterior boundaries as thermally insulated.

The pressure in the cavity after a PD event can be calculated using the ideal gas law:

$$p_{tPD+\Delta tPD}=\left(\frac{T_{tPD+\Delta tPD}}{T_{tPD}}\right)p_{tPD},$$

(10)

where $p_{tPD}$ and $p_{tPD+\Delta tPD}$ are, respectively, the pressure in the cavity before and after the PD event. Similarly, $T_{tPD}$ and $T_{tPD+\Delta tPD}$ are, respectively, the temperature in the cavity before and after the PD event.

3. Case of Study and Simulation Results

The simulation model was implemented as follows: the analytical model, Expressions (1)–(7) were implemented using a code in MATLAB. On the other hand, the temperature model in Equation (8), was solved using the PDE toolbox of MATLAB which implements a finite element analysis method and uses, in turn, an internal solver based on the variable order method [28]. The absolute and relative tolerances of the solver are, respectively, $1\times {10}^{-6}$ and $1\times {10}^{-3}$. Finally, the pressure was calculated for each time step using Equation (10).

Simulations were implemented for a considerable number of cycles (100 at phases A through C, 200 at phase E and 300 at phase C) in order to consider the stochastic variations in the model. A time step of $1\times {10}^{-5}$ s was used for executing simulations when there are not PD acting (inception and electron generation rate criteria are not satisfied) and $1\times {10}^{-9}$ s during PD events.

The same case used in [23] is considered in this study. However, in order to analyze the temperature and pressure behavior after PD activity and its effects on PD behavior, Equations (1)–(3) and (5), are dynamically calculated along the simulation time. The geometry consists of a linear, homogenous and isotropic dielectric bulk of epoxy resin of 3.5 mm thickness between two parallel plates, and a spherical void, filled with air, of 2.5 mm diameter is included in the middle of the dielectric bulk. A 19.25 kV, 50 Hz sinusoidal voltage source is applied to the upper electrode while the lower is grounded.

Table 1. Parameters of case of study. Thermal parameters taken from [20].

Parameter	Description	Value
$a$	Cavity radius	1.25 mm
$h_m$	Homogeneous material thickness	3.5 mm
${\epsilon }_{r\_mat}$	Relative Permittivity of dielectric material	4
${\epsilon }_{r\_cav}$	Relative Permittivity of cavity	1
$C_{p\_mat}$	Material heat capacity	1179 J·kg−1·K−1
$C_{p\_air}$	Air heat capacity	1005 J·kg−1·K−1
$k_{T\_mat}$	Material thermal conductivity	0.35 W·m−1·K−1
$k_{T\_air}$	Air thermal conductivity	0.0257 W·m−1·K−1
${\rho }_{m\_mat}$	Material volume density	561 kg·m−3
${\rho }_{m\_air}$	Air volume density	1.205 kg·m−3

summarizes the parameters of media related to the case of study.

In order to make comparisons with measurements presented in [23], it is considered that cavity surface conductivity is 0 S·m⁻¹.

Five different consecutive ageing phases were considered, namely: phase A, with a duration up to 10 min; phase B, lasts about 20 to 50 h; phase C, lasts about 100 to 200 h; phase D, lasts about 1200 h and phase E, lasts 50 h. Measurements and simulations using the analytical model were reported in [23] for each of the above ageing phases. Figure 2 shows the phase-resolved partial discharge (PRPD) pattern and q-φ-n diagrams obtained for the case of study at the five considered ageing phases.

Table 2. Parameters for the PD model simulation at each ageing phase, taken from [23], and simulation results summary for PD charge magnitude.

Parameter	Phase A	Phase B	Phase C	Phase D	Phase E
${\Phi }_{dt}$ (eV)	1	1.1	0.98	1.45	1.33
$\tau$ (s)	0.002	0.002	0.002	1000	1000
$p_0$ (kPa)	65	65	6	2	2
$T_0$ (K)	300	300	300	300	300
${\xi }_{+}/{\xi }_{-}$	50	50	50	50	1
$\left|q_{\max }\right|$ (nC)	0.809	1.585	0.171	5.519	3.753
$\left|q_{\min }\right|$ (nC)	0.800	0.800	0.162	0.092	0.085
$\left|q_{mean}\right|$ (nC)	0.805	0.824	0.166	3.438	1.300
NHW simulated	6	5.6	33	0.03	2.2
NHW measured [23]	5–7	~5.8	~30	0.03–0.06	~3
NHW error (%)	14.286	3.448	10	0	26.667
Number of simulated cycles	100	100	100	300	200
$E_{incm}$ (kV·mm−1)	2.642	2.641	0.468	0.235	0.235
$E_{extm}$ (kV·mm-1)	0.592	0.592	0.055	0.018	0.018

exhibits the parameters of the stochastic model and initial conditions used for simulating each ageing phase, as well as a summary of the main results for the PD electric charge magnitude.

In Table 2, $E_{incm}$ and $E_{extm}$, respectively, are the average inception and extinction-field strength magnitudes along the simulation time, N_HW is the number of PD per half cycle, and $p_0$ and $T_0$ are, respectively, the initial pressure and temperature of gas inside the cavity.

Simulation results show reasonable agreement with measured values reported in [23]. When compared with measured values, the simulated number of PD per cycle (N_HW) present a maximal difference of 26.667% at phase E (N_HW measured ~3 and simulated 2.2). However, at this phase it must be considered that the pressure behavior with time has a complex dependence on temperature and applied voltage that the ideal gas law cannot represent completely [29].

Figure 3 shows the temperature and pressure simulation results for the case of study during the first five periods of the AC-applied voltage at ageing phases A, B and E.

In Figure 3, Tmcav is the average temperature in the cavity volume and Tscav is the average temperature at the void surface. Taking into account that maximal changes in pressure and temperature at phases C and D are below 1 K and 0.01 kPa, they were not plotted.

Table 3. Temperature and pressure, summary of simulation results at each ageing phase.

Variable	Phase A	Phase B	Phase C	Phase D	Phase E
Max $T_{mcav}$ (K)	302.7	305.0	300.3	347.4	321.9
Min $T_{mcav}$ (K)	300.0	300.0	300.0	300.0	300.0
Mean $T_{mcav}$ (K)	301.0	301.0	300.2	300.1	300.9
Max $T_{scav}$ (K)	301.5	302.9	300.2	327.2	312.6
Min $T_{scav}$ (K)	300.0	300.0	300.0	300.0	300.0
Mean $T_{scav}$ (K)	300.6	300.6	300.1	300.0	300.5
Max $p$ (kPa)	65.6	66.1	6.0	2.3	2.2
Min $p$ (kPa)	65.0	65.0	6.0	2.0	2.0
Mean $p$ (kPa)	65.2	65.2	6.0	2.0	2.0

shows the summary of simulation results for the five ageing phases considered and simulation periods in Table 2.

4. Discussion

At phase A, the PRPD pattern exhibits a “bar-like” structure in which the PD charge magnitude is almost constant with an average magnitude very close to the minimum value. The measured minimum magnitude and number of PD per half cycle were, respectively, 0.9 nC and about 5 to 7 [23]. Simulated values in the second column of Table 2 are in good agreement with measured values. On the other hand, in Figure 3a, it can be seen that after each PD event, the temperature in both the gas in cavity and the void surface, increases. However, the average increase in temperature along the simulated periods is just 1.0 K. A similar behavior is found for pressure, shown in Figure 3b, where the increase in average pressure is 0.2 kPa.

At phase B, Figure 2c, the PRPD pattern exhibits a “bar-like” plus a “rabbit ear-like” structure. The “bar-like” structure is close to around the minimum PD charge magnitude, while the tip of the “rabbit ear-like” structure tends to increase. For the number of simulated cycles, the maximum PD charge magnitude calculated was 1.6 nC. In order to reproduce the experimental measurements at this phase, the detrapping work function is increased to 1.1 eV. This produces an increase in the time lag, and this in turn, a reduction in the number of PD per half cycle and an increase in the PD charge magnitude. The increased time lag also means that more energy stored is converted into heat at each PD event, which can be seen at Table 3. However, the average temperature in cavity gas and cavity surface is reduced in comparison to phase A, because the time for cooling after a PD event is increased if the number of PD per half cycle is reduced. A similar behavior is found in pressure; the increment of maximum pressure magnitude is negligible.

In order to reproduce experimental measurements at ageing phases C to E, the initial pressure has to be decreased for considering the changes in the gas consumption and gaseous by-products generation rates [30]. After about 30 min, in ageing phase B, the pressure in the cavity starts to decrease continuously [18] with a decreasing rate that is voltage and temperature dependent [29].

At phase C, Figure 2e, the PRPD pattern exhibits a “bar-like” structure similar to phase A. However, the magnitude of the PD charge decreases in comparison to phase A by 20.2%, and the number of PD per half cycle increases up to 550.4% compared with phase A. On the other hand, the q-φ-n diagram exhibits an almost continuous distribution, which means an almost constant PD condition in the cavity. The latter can be seen also in Figure 3e,f where the average temperature is continuously changing. The initial cavity gas pressure is reduced to 6 kPa, which involves a decreasing in the inception and extinction field magnitudes; this produces a decrement in the time lag and a reduction in the PD charge magnitude. In spite of almost continuous PD activity inside the cavity, the temperature and pressure in the cavity remains almost constant, and the average increment in temperature and pressure is just 0.1% when it is compared with room and initial conditions.

At phase D, Figure 2g, the PRPD pattern exhibits an irregular behavior with marked intermittency in the PD time lag. The PD charge magnitude exhibits a great increase with regard to phase A of 581.7% while the number of PD per half cycle is just 0.03. The measured number of PD per cycle is in the range 0.02 to 0.06, which means that simulation results are in reasonable agreement with values measured in [23]. In order to reproduce experimental measurements, the detrapping work function and the effective detrapping time constant are increased while the initial pressure is set to 2 kPa. With these changes, the time lag is hardly increased. Despite the reduction of inception magnitude, the reduction of the electron generation rate dominates the stochastic PD behavior. The average temperature and pressure exhibit a light increase of 0.03%. Despite the increase of energy transformed into heat, the maximum temperature increases in 47.4 K compared with room temperature, and the larger time between consecutive PD eases a rapid cooling.

At phase E, Figure 2i, the PD behavior is continuous and the PD charge magnitude seems to be uniformly distributed between two sinusoidal envelopes in the PRPD pattern. The magnitude of the PD charge is reduced in comparison to phase D by 32%. However, the number of PD per half cycle is increased up to 2.2, which is in reasonable agreement with the experimental value of about 3 [23]. Pressure and temperature exhibit similar behavior and magnitudes as for phase D, Figure 3i,j, with slight differences in average temperature and pressure, as shown in Table 3.

Despite the agreement with experimental measurements, the PD behavior at advanced ageing conditions, such as phase E, cannot be explained only considering the changes of the gas filling the cavity. The physical and chemical interactions at the solid—gas interface must be considered.

In this study, the effect of cavity conductivity on PD behavior and its variation along the ageing process is neglected in order to analyze the effects of temperature and pressure variation on the stochastic behavior; however, cavity conductivity plays a fundamental role at PD activity and is related to surface conditions, and for that reason, in future works it will be considered.

Simulations were carried out using a MATLAB application implemented by authors named PDSym1S developed for these purposes.

5. Conclusions

Temperature and pressure changes inside an air-filled cavity within a homogeneous dielectric material after PD events were calculated using an analytical-FEA model. Simulation results exhibit reasonable agreement with experimental measurements. Pressure in the gas filling the cavity is one of the main parameters affecting the PD behavior at ageing phases. Taking into account that the highest temperature increase along the simulated phases was only about 15%, temperature effect on surface degradation can be neglected in comparison to chemical reactions and particle impact, which is in line with other experimental studies. Future work will include the consideration of different electrodes geometries and unevenness of electrodes surfaces.