Partial Discharge Modelling Based on a
Cylindrical Model in Solid Dielectrics
George Chen and Fauzan Baharudin
School of Electronics and Computer Science
University of Southampton, United Kingdom
Email: gc@ecs.soton.ac.uk
Abstract--Partial discharges (PDs) are harmful to electrical
insulation as they may degrade the material. Consequently, a lot
of efforts have been made to develop suitable systems for PD
measurements and monitoring. Due to the complexity of PDs and
many factors that can influence PD characteristics, the research
on PDs is still very active and several models have been proposed
over the years to understand the behaviours of PDs.
In the present paper, we have investigated the frequency
characteristics of a cylindrical void in solid dielectric based on
the model proposed by Forssen. Using similar program
environments of MATLAB and COMSOL, our algorithm
performs better in computation time and gives the similar results
compared with the experiments. More importantly, a new
interface has been developed so users are able to change
parameters. This new feature allows the user to examine effects
of various factors such as conductivity of the cavity surface, the
applied voltage and statistical time lag on PD characteristics.
Index Term--Partial discharge, simulation, solid dielectrics,
frequency effect.
I. I
NTRODUCTION
a
d
rtial discharges (PD) are localised electric discharges that
o not bridge the whole distance between electrodes;
indicating the presence of cavity and defects in the insulation
such as high voltage power cables and power generators. The
measurement and investigation of partial discharges in cavity
to diagnose the condition of the insulation system has been of
interest since the early 1940’s. PDs are generally divided into
three different groups because of their different origins: (1)
corona discharges, (2) internal discharges and (3) surface
discharges [1].
Traditionally, PDs are only measured and examined at a
single frequency of the applied voltage. However, it is now
possible to measure partial discharges at variable frequencies
of the applied voltage using a rather new technique called
Variable Frequency – Phase Resolved PD Analysis (VF-
PRPDA) method [2]. The new technique provides researchers
more information and data about the state of the insulation
system as compared to a single frequency.
In the present paper, we have investigated the frequency
characteristics of a cylindrical void in solid dielectric based on
the model proposed by Forssen [3]. Using similar program
environments of MATLAB and COMSOL, a new interface
has been developed so users are able to change parameters.
This new feature allows the user to examine effects of various
factors such as conductivity of the cavity surface, the applied
voltage and statistical time lag on PD characteristics.
II. PD
MODEL
In laboratory investigation a sandwich model is often used
as shown in Figure 1. For comparison we use the same
arrangement as in [3]. Three polycarbonate plates each of 1
mm thick and 14 mm in diameter are placed on top of each
other with the middle plate being drilled at the centre to
resemble a void with 10mm of diameter. A 10 kV sinusoidal
voltage was initially applied between the electrodes to
investigate the influence of frequency, ranging from 0.01Hz to
100Hz. To prevent corona discharges the sample and
electrodes are typically casted with epoxy resin.
Po l ycar b o n at e
plates
Epo xy r esi n
Ca v i t y
Brass elect rode
Figure 1 Illustration of the PD model simulated.
The basic governing equations are as below:
f
D
ρ
∇
⋅= （1）
0
f
f
J
t
ρ
∂
∇
⋅+ =
∂
（2）
with D as the electric displacement field, ρ
f
as the free charge
density, and J
f
as the free current density. But for non-
dispersive, linear isotropic dielectrics exposed to slowly
varying fields, equations 1 and 2 can be combined to obtain
equation 3; where σ equals to the electric conductivity, V
equals to the applied electric potential and ε as the permittivity.
P
978-1-4244-1622-6/08/$25.00 ©2008 IEEE
2008 International Conference on Condition Monitoring and Diagnosis, Beijing, China, April 21-24, 2008
(()
V
V
t
σε
∂
∇⋅ − ∇ − ∇ =
∂
)0
（3）
Partial discharges will occur if and only if the insulation
systems satisfy two necessary conditions: (1) There must be
free electrons available at the surfaces of the void to start an
electron avalanche process and (2) The electric field must be
high enough, corresponding to the associated voltage called
the Inception Voltage, U
inc
.
A. Statistical Time Lag
The period between U
inc
and the starting point of the
electron avalanche due to the delay in getting free electrons is
called Statistical Time Lag, τ
stat
. The discharges will then
continue until the electric field in the cavity becomes too low
such that it reaches a specific voltage called the Extinction
Voltage, U
ext
.
If τ
stat
is negligible compared to the period of the applied
voltage, T
v
, it will not be included in the simulation to speed
up processing. Up until now however, there is no clear
distinction between significant and negligible τ
stat
value and
this imposes problems when modelling PD activity. Therefore,
the following criteria have been introduced in this simulation
to determine the borderline between significant and negligible
value for PD analysis as below:
10000
10000
v
stat
v
stat
T
Negligible
T
Significant
τ
τ
≥
<
（4）
B. Frequency dependent nature of PDs
There are also two time constants that affect and
influence the PD activity: (1) Charge redistribution time on
the cavity surface, τ
cavity
. Determined from the conductivity of
the cavity surface, σ
surf
, and the geometry of the void; e.g. a
higher σ
surf
and/or a smaller void will results in a shorter τ
cavity
.
(2) The charge redistribution time in the surrounding
insulation material, τ
material
. Determined from the conductivity
and permittivity of the bulk (surrounding) insulation material
as well as its geometry; e.g. an increase in conductivity and/or
a decrease in permittivity will results in a shorter τ
material
.
PD activity is dependent on the frequency of applied
voltage only when:
[/]
[]
cative material v
stat v
and or or T and or
or T
ττ
τ
<∋
>∋
/;
（5）
where ∋ means “within the same range as”. Both τ
cavity
and
τ
material
value can be associated with the ‘Screening Effect’ and
‘Blocking Effect’ depending on mutual relationship between
them.
PD activity in the cavity was modelled as an increase in the
associated channel conductivity value and they were modelled
to behave as equation 6 below; where σ
0
is the initial
channel’s electrical conductivity value before discharging, I is
the current across the channel, U is the electric potential
difference across the channel and I
crit
is the critical current
value to start an electron avalanche [3]. For numerical reasons
the conductivity needs to be limited to a maximum value of
σ
max
(i.e. 1x10
-4
) in order to achieve numerical convergence
for the final solution and calculation.
arg
int
0
0
vI
during disch e
vI
crit
e
otherwise
σ
σ
σ
⎛⎞
⎜+ ⎟
⎜⎟
⎝⎠
⎧
⎪
=
⎨
⎪
⎩
（6）
{
arg
0
ss
J
dS VdS during disch e
otherwise
I
σ
∫∫ =∫∫ ∇
=
（7）
Equation 7 above is the equation for the current through the
discharging channel [3]. Notice that I and σ in both of the
equations above are dependent on each other. The only
difference is that they are of a different time step; i.e. σ takes
the value of I from the previous time step.
At the frequency where τ
stat
is significant in comparison to
the T
v
, equation 8 below needs to be used in addition to the
U
inc
limit alone to initiate PD [3]:
0
0
U
U
inc
e
inc
inc
Ie
dt U U
N
UU
P
>
<
⎧
⎪
⎪
=
⎨
⎪
⎪
⎩
(8)
where P is the probability of PD occurring, I
e0
is the initial
electron injection intensity and N is the number of channels
inside the cavity. The value of dt needs to be calculated using
equation 9:
0
int
360
time ste
pp
in
g
erval durin
g
no PDactivit
y
dt
frequency
=
×
（9）
so that P can be compared with a random number R which lies
within the range of 0 < R < 1 such that [3]:
（10）
PRPDoccurs
P R no PD occurs
>
<
For the case of having more than one channel that exceeds
U
inc
, only one channel will be chosen randomly to discharge
by MATLAB.
C. COMSOL and MATLAb in Brevity
The COMSOL software used in this project based entirely
on the finite element method (FEM) structure analysis. It
solves various physical problems by finding the approximate
solution of partial differential equations (PDEs) and then
renders them into equivalent ordinary differential equations
(ODEs) to be solved using standard techniques such as the
Newton-Raphson iteration method [4]. The full equation used
by COMSOL in solving the electromagnetics module is as
below:
2
2
()
aa
VV
e d c V aV aV V f
tt
γβ
∂∂
++∇⋅−∇−+++∇=
∂∂
(11)
When the entire zero valued variable from Table are
eliminated, the resulting equation is simplified and becomes:
0
()
r
djw V0
σ
εε
−∇⋅ + ∇ = (12)
where d is the thickness of the sample, ω is the frequency in
radians and ε
0
as the permittivity of vacuum and ε
r
as the
relative permittivity.
Table I
Variable definitions for equation 11
a = Absorption coefficient = 0
f = Source term = 0
e
a
= Mass coefficient = 0
d
a
= Damping coefficient = 0
α = Conservative flux convection
coefficient
= 0|0 (matrix)
β = Convection coefficient = 0|0 (matrix)
γ = Conservative flux source term = 0|0 (matrix)
c = Diffusion coefficient = d.(σ+jωε
0
ε
r
)
MATLAB on the other hand combines both numerical
computing environment and programming language. Using
matrix manipulations, MATLAB allows the user to easily plot
functions or data, implement algorithms, and create graphical
user interfaces (GUI). With seamless interoperability between
COMSOL and MATLAB, they have been widely used by
researchers and designers to model, simulate and solve various
physical problems in the engineering field.
The model was drawn under the 2D axial symmetry space
dimension using the AC/DC module within the COMSOL. To
add some random behaviour of the PD activity, the void was
divided into 5 different channels; each having equal volume of
space. This was done to illustrate that the discharges may
occur at only some parts of the cavity. Cavity surface was
added in the model by adding layers of 0.1mm thick
surrounding the void.
D.
COMSOL with MATLAB Interface
Iterative elements were introduced in the model using
MATLAB. The coding was structured initially using low
frequency value of the applied voltage and then followed by a
higher frequency; mainly because of the τ
stat
factor. COMSOL
provides an excellent interoperability with MATLAB.
Figure 2 Interactive menu at the beginning of the simulation with default
parameter values
Initial analysis revealed that running 20 cycles is sufficient
to achieve a respectable result and feasible simulation time. To
increase simulation speed further, the variables were pre-
allocated initially using the ‘zeros’ function wherever possible;
especially within a loop (the coding itself contains more than
20 variable names). In addition, the whole coding structure
was made on the basis of having minimum conditional
syntaxes and iterations where possible.
To investigate the effect of the amplitude the applied
voltage the voltage was varied from 8kV to 12kV in 2kV steps
to see any effects towards the PD activity inside the channels.
The σ
surf
was also varied from 1x10
-10
up to 1x10
-20
S/m to see
any variations in the PD activity. Most of the additional
investigations were simulated once at the low frequency
region and once again at the high frequency. With the help of
interactive menu at the start of the simulation as shown in
Figure 2, the variables used for the additional investigations
can be keyed in easily and quickly.
III. R
ESULTS AND DISCUSSIONS
A. Average PD per Cycle
The 0.1 Hz frequency took the longest time of
approximately 16 hours to finish simulating the model
whereas the 100 Hz took the shortest duration of 11 hours. It
might be safe to assume generally that the time it takes to
simulate the model for each frequency is inversely
proportionate with the frequency of the applied voltage.
However, it does not necessarily means that the PD activity
will also follow this trend and this is reflected in the 0.01 Hz
data where the duration is in fact shorter than the next higher
frequency which is the 0.1 Hz. To explain this further, refer to
Figure 3.
At 0.01Hz, its average PD per cycle is in fact lower than
the average PD per cycle of the 0.1Hz frequency. This means,
the total simulation time per frequency is behaving
proportionally to the number of discharges per cycle; and this
seems to be logical because for every PD occurrences,
MATLAB will switch into the PD activity program loop
which uses 0.02° as the time stepping interval. For example,
instead of stepping from 45° to 48° (default 3° interval),
MATLAB will sweep in 0.02° time step from 45° up to an
angle at which the electric potential across the specific
discharging channel reaches U
ext
(usually after 3 times of
0.02° time step i.e. 45.06°). Only then MATLAB will revert
back to the 3° time steps and repeat the simulation again
starting at 48° phase angle. In other words, COMSOL will be
called at least three times more to account for the additional
time stepping interval.
At frequencies above 0.1Hz (i.e. where the τ
stat
is
significant), the program will simulate τ
stat
influence to the PD
activity by including equation 8, 9 and 10 in choosing only
one to discharge for the next time step. Clearly, the 1Hz
frequency is affected by the inclusion of those additional
equations and hence extra time spent. Although the additional
time might not be much, but over 20 cycles, the time delay
sums up to cost 1Hz frequency additional 30 minutes; longer
than 0.01Hz even though 0.01Hz has more PD.
Figure 3 Average PD per cycle for various frequencies: lab experiment,
existing PD simulation model and present simulation.
Referring to the trend line in Figure 3, it can be seen that
close agreements have been reached between the simulated
model, Forssén’s model and the actual practical experiment
[3]. This ultimately proves that the simulation model produced
and programming code generated for this project is correct
and acceptable as a whole to be used for any other PD analysis
using different parameter values.
B. τ
cavity
and τ
material
Influences
Altering σ
surf
value is one way of investigating the
relationship between τ
cavity
and τ
material
and their influences
towards PD activity because τ
cavity
behaves inversely
proportional to σ
surf
. Simulations were done for σ
surf
using
values of 1x10
-10
S/m, 1x10
-15
S/m (default) and 1x10
-20
S/m at
0.1Hz and 50Hz of the applied frequency. The insulation
conductivity value, σ
ins
, was fixed at 1x10
-15
S/m throughout
all simulations. Table 2 summarises the simulated results for
0.1Hz and 50Hz.
Table II
Average number of PD per cycle using various cavity surface conductivity
values; simulated at 0.1Hz and 50Hz frequency.
Average PD per cycle Applied
frequency (Hz)
1x10
-20
S/m 1x10
-15
S/m
(default)
1x10
-10
S/m
0.1
10.05 9.9 9.9
50
4.4 3.8 3.6
There will be no difference at all in terms of physical
characteristics when the cavity surface and the surrounding
insulation material are defined using the default values for
their electrical conductivity and relative permittivity (refer to
Figure 2); i.e. (τ
cavity
= τ
material
). However, increasing σ
surf
to
1x10
-10
S/m (to simulate ageing) will effectively decrease the
value of τ
cavity
; turning the relationship into (τ
cavity
< τ
material
).
This causes the ‘Screening Effect’ to occur which effectively
reduces the average number of PD per cycle. The opposite is
true for reducing σ
surf
; causing the ‘Blocking Effect’ to occur
as mentioned in the Background Reading section. Table
confirms these relationships as it suggested that the average
number of PD per cycle will decrease as the σ
surf
increases.
C. Altering the Amplitude of the Applied Voltage
At low amplitude e.g. 8000V, the average time it took for
any of the channels to stay above the U
inc
limit will be much
less than the one experiencing high amplitude of the applied
electric potentials e.g. 12000V. This directly decreases the
chances for any of the channels to discharge and thus the
average PD per cycle will reduce. Conversely is true for the
high amplitude of the applied voltage (refer to Table 3).
Table III
The effect of altering the amplitude of the applied voltage; simulated at 50Hz
Amplitude of the applied voltage (V) at 50Hz
80000 10000 12000
Average PD per
cycle
1.75 3.8 5.45
IV. CONCLUSIONS
COMSOL and MATLAB work seamlessly with each other,
making it possible to simulate the complex and random nature
of PD activities. The simulation platform established in the
present work is proven to be ideal for simulating PD activity
in an otherwise homogeneous dielectric material exposed to
high electric fields. This simulation platform is better than the
existing model in the following two aspects:
In the present simulation platform, the value dt changes as
the frequency changes, making it versatile for all simulation
work and more efficient.
The interactive menu introduced is useful for further
investigations using various material properties or PD
parameter values.
Specifically for the PD activity, it can be summarised that
the τ
stat
value affects the number of average PD per cycle, at
the frequency above 0.1Hz i.e. where the τ
stat
is found be
significant to the T
v
, determined using equation 4.
Varying the σ
surf
value will in effect alter the value of τ
cavity
and in turn observes the ‘Screening Effect’ or the ‘Blocking
Effect’. Screening occurs when (τ
cavity
< τ
material
) and blocking
occurs when (τ
cavity
> τ
material
).
Lastly, lowering the amplitude of the applied voltage will
reduce the average number of PD per cycle due to the reduced
area above the U
inc
line and vice-versa.
V. REFERENCES
[1] F. H. Kreuger, Discharge Detection in High Voltage Equipment,
Butterworths, London, 1989.
[2] Hans Edin, Partial Discharges Studied with Variable Frequency of the
Applied Voltage. Sept 14, Stockholm, Sweden : Kungl Tekniska
Högskolan, Vols. TRITA-EEK-0102. ISSN 1100-1593, 2001.
[3] H. Edin and C. Forssen. Modelling of a discharging cavity in a dielectric
material exposed to high electric fields. KTH IR-EE-EEK. 2005, 2005:003.
[4] J. H. Mathews and K. D. Fink, Numerical Methods: Using Matlab, Fourth
Edition, Prentice- Hall Pub. Inc, 2004.
