0 Bewertungen0% fanden dieses Dokument nützlich (0 Abstimmungen)

7 Ansichten137 SeitenThis shows the modeling and simulation of varistors.

© © All Rights Reserved

PDF, TXT oder online auf Scribd lesen

This shows the modeling and simulation of varistors.

© All Rights Reserved

Als PDF, TXT **herunterladen** oder online auf Scribd lesen

0 Bewertungen0% fanden dieses Dokument nützlich (0 Abstimmungen)

7 Ansichten137 SeitenThis shows the modeling and simulation of varistors.

© All Rights Reserved

Als PDF, TXT **herunterladen** oder online auf Scribd lesen

Sie sind auf Seite 1von 137

of Varistors

Modellierung und Simulation von Varistoren

Vom Fachbereich Elektrotechnik und Informationstechnik der Technischen

Universität Darmstadt zur Erlangung des akademischen Grades

Doktor-Ingenieur (Dr.-Ing.) genehmigte Dissertation von Dipl.-Wirtsch.-Ing.

Frank Stefan Denz aus Frankfurt am Main

Tag der Einreichung: 22.04.2014, Tag der Prüfung: 08.09.2014

Darmstadt 2014 — D 17

2. Gutachten: Prof. Dr.-Ing. V. Hinrichsen

Modeling and Simulation of Varistors

Modellierung und Simulation von Varistoren

am Main

2. Gutachten: Prof. Dr.-Ing. V. Hinrichsen

Tag der Prüfung: 08.09.2014

Darmstadt — D 17

URN: urn:nbn:de:tuda-tuprints-41971

URL: http://tuprints.ulb.tu-darmstadt.de/4197

E-Publishing-Service der TU Darmstadt

http://tuprints.ulb.tu-darmstadt.de

tuprints@ulb.tu-darmstadt.de

Namensnennung – Keine kommerzielle Nutzung – Keine Bearbeitung 2.0 Deutsch-

land

http://creativecommons.org/licenses/by-nc-nd/2.0/de/

Erklärung laut §9 PromO

Ich versichere hiermit, dass ich die vorliegende Dissertation allein und

nur unter Verwendung der angegebenen Literatur verfasst habe. Die

Arbeit hat bisher noch nicht zu Prüfungszwecken gedient.

(Frank Denz)

1

Kurzfassung

Diese Arbeit befasst sich mit verschiedenen Problemen im Zusammenhang mit Va-

ristoren und Mikrovaristoren, welche dank ihrer außergewöhnlichen, nichtlinearen

elektrischen Leitfähigkeit zur Unterdrückung transienter Überspannungen eingesetzt

werden. Diese Arbeit ist vor allem dadurch motiviert, dass man zum einen das Ver-

halten von Überspannungsableitern im Hochspannungsbereich simulieren möchte,

zum anderen dasjenige von Mikrovaristoren für zukünftige Einsatzmöglichkeiten

bei der nichtlinearen resistiven Feldsteuerung.

Die Untersuchung der Überspannungsableiter erfordert die numerische Berech-

nung wechselseitig gekoppelter elektrischer und thermischer Felder, wobei die

Hauptschwierigkeit in der extremen Nichtlinearität des elektrischen Teilproblems

zu finden ist. Zu diesem Zweck wird die Elektroquasistatik-Gleichung mittels der

Finite-Elemente-Methode im Zeitbereich gelöst. Auf die Berechnung des thermisch-

stationären Zustandes eines Überspannungsableiters und die Untersuchung eines

Enveloppengleichungsmodells zur Simulation des Erwärmungs- und Kühlverhaltens

von Ableitern wird näher eingegangen.

Diese Simulationen sind abhängig von hinreichend genauen Modellen zur Be-

schreibung der Materialeigenschaften. Die Schätzung der nichtlinearen Leitfähigkeit

und Permittivität von Varistormaterialien ist ein essenzieller Bestandteil dieser

Arbeit.

Des Weiteren werden nichtlineare Kapazitäts- und Leitfähigkeitsmatrizen einge-

führt. Der hier vorgestellte Ansatz beruht auf einem Ersatzschaltungsmodell, dessen

Parameter mithilfe von Feldsimulationsergebnissen bestimmt werden.

Abstract

This thesis treats various problems that arise in the context of varistors and

microvaristors, which are are used for the suppression of transient overvoltages, due

to their extraordinary nonlinear electrical conductivity. The present work is mainly

motivated by the desire to simulate the behavior of high-voltage surge arresters

used for lightning protection on the one hand and of microvaristors as materials for

future applications in nonlinear resistive stress control on the other hand.

The analysis of surge arresters requires the numerical calculation of mutually-

dependent electric and thermal fields, whereby the principal difficulty resides in

the extreme nonlinearity of the electric problem. For this purpose, the electro-

quasistatics equation is solved in time domain by means of the finite-element

method. The calculation of the thermally stationary state of a surge arrester and the

3

evaluation of an envelope equation model for simulating the heating and cooling

behavior of arresters are discussed in more detail.

These simulations depend on sufficiently accurate models that describe the mate-

rial properties. The estimation of nonlinear conductivity and permittivity of varistor

materials is an inherent part of this thesis.

Furthermore, nonlinear capacitance and conductance matrices are introduced.

The presented approach is based on an equivalent circuit model. Its parameters are

determined from field-simulation results.

4

Contents

1 Introduction 7

3 Physical Background 17

3.1 Maxwell’s Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.1.1 Field Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.1.2 Constitutive Equations . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1.3 Electro-Quasistatic Approximation . . . . . . . . . . . . . . . . . 21

3.2 Heat Conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.2.1 Heat Conduction Equation . . . . . . . . . . . . . . . . . . . . . . 25

3.2.2 Convection and Radiation Conditions . . . . . . . . . . . . . . . 26

3.3 Metal-Oxide Varistors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.3.1 Microstructure, Manufacturing and Destruction Mechanisms . 29

3.3.2 Constitutive Relations . . . . . . . . . . . . . . . . . . . . . . . . . 32

4.1 Introduction to the Finite Element Method . . . . . . . . . . . . . . . . 35

4.2 Application to Electro-Quasistatic Problems . . . . . . . . . . . . . . . . 38

4.2.1 Discrete Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4.2.2 Solution in Time Domain . . . . . . . . . . . . . . . . . . . . . . . 40

4.2.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.3 Application to Heat-Conduction Problems . . . . . . . . . . . . . . . . . 45

4.3.1 Discrete Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.3.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.1 Nonlinear Capacitances . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.2 Measurement Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

5.3 Conventional Characterization Method . . . . . . . . . . . . . . . . . . . 52

5.4 Least-Squares Based Method . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.4.1 Description of the Method . . . . . . . . . . . . . . . . . . . . . . 54

5.4.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

5

5.4.3 Characterization of a Microvaristor-Filled Silicone-Rubber Ma-

terial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.5 Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.1 Computation of Stationary Operational Condition of Surge Arresters 76

6.2 Envelope Model for Transitory Processes in Surge Arresters . . . . . . 84

7.1 Capacitance and Conductance Matrices . . . . . . . . . . . . . . . . . . 99

7.2 Nonlinear Matrices and Equivalent Circuits . . . . . . . . . . . . . . . . 102

7.3 Extraction of Nonlinear Circuit Parameters . . . . . . . . . . . . . . . . 103

7.4 Nonlinear Equivalent Circuit of Multi-Conductor Cable . . . . . . . . . 106

Appendix 115

Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Bibliography 125

Acknowledgements 135

6

1 Introduction

In the last few decades scientists and engineers have increasingly turned towards

computational methods to find solutions for their problems. This has been motivated

by the continuous advance in computer technology. Over more than five decades the

power of computers has tended to increase exponentially [68, 77]. Simultaneously,

the cost of computer systems has come down to the extent that every mobile phone

or washing machine can perform calculations, which few could have imagined a few

years ago. When Argyris laid the foundations of the finite element method during

World War II, he initially disposed of an electromechanical calculator to solve a

structural design problem of an airplane with 64 linear equations [97]. Today it is

possible to solve problems with hundreds of millions or more unknowns [80].

In the course of time, the numerical simulation of electromagnetic fields has

become an important tool in many areas of electrical engineering from high fre-

quency to low frequency, e.g., the design of machines and other components which

are essential for generation, transport and utilization of electrical energy. In many

cases, it is possible to reduce the number of experimental measurements or even to

completely replace measurements by simulations. Generally, computer simulations

require less time and money than the construction of new prototypes leading to

faster and less expensive product development. Simulations may also offer insights,

which can not be obtained from measurements. This may be the case, if the physical

quantity, in which one is interested, is inaccessible from the outside or if the mea-

surement device significantly affects the outcome of the measurement. Simulations

are also nondestructive and reproducible, i.e., not influenced by uncontrolled factors.

Finally, some experiments are avoided for ethical reasons (e.g., biological effects of

radiation) or simply impossible.

However, there are still many problems which require a more traditional ap-

proach. The number of problems, for which perfectly appropriate methods exist,

but which are simply ’too large’ to be treated numerically, diminishes continuously

as memory size and speed of computers increase. Yet, many problems remain

outside of the scope of numerical methods either because the physical properties

needed for a model are unknown or because the underlying physics can not be

handled by readily available software. In these cases, trial & error coupled with

experience and an intuitive understanding remains of greatest importance. Shifting

the frontier between the experimental and simulation approaches by extending the

7

range of problems which can be treated numerically is one of the principal goals in

computational electrodynamics.

This work is dedicated to a particular class of problems, the modeling and sim-

ulation of varistor-related problems. The initial objective was the simulation of

high-voltage surge arrester behavior during lightning-surge impulses. The main

functional component of such surge arresters are metal-oxide varistors (MOVs).

Their behavior is not independent from temperature. During a lightning surge

the thermal energy increases significantly. Therefore, the choice of a multiphysics

approach combining electromagnetics, more precisely electro-quasistatics, with heat

conduction is indispensable. However, the main challenge is the extreme nonlinear-

ity of the varistor conductivity with respect to field strength. This nonlinearity leads

to time steps, which are so small that the simulation of a lightning surge of several

microseconds may take weeks, even when employing a computer cluster. Due to the

high costs in computational power and time to obtain results of limited practical

value, attention passed to other varistor-related problems, mainly connected with

research group FOR 575 of the Deutsche Forschungsgemeinschaft (DFG). Ultimately,

the author’s work on the transient multiphysics simulation of lightning surge pulses

is completely omitted from this document [16, 17].

The topics of research group FOR 575 include the research on means to filter

traveling waves emanating from insulated-gate bipolar transistor (IGBT) inverters.

One approach is the embedding of microvaristor materials into the insulation of

wires to achieve a nonlinear resistive stress grading. As microvaristors are not

fundamentally different from normal varistors, they can be simulated in the same

way as normal varistor materials. However, it became apparent that the knowledge

of their material properties was insufficient and that, prior to any simulation, one had

to find a way to improve the quality of the material models. For that reason, a novel

approach has been developed for the characterization of the nonlinear electrical

properties, which is described in Ch. 5 along with the resulting observations for a

specific microvaristor material.

While the transient simulation of surge arrester behavior during short current

pulses is omitted, two more specialized kinds of surge arrester simulation are

presented in Ch. 6. The first is concerned with the coupling of a transient electrical

simulation with static heat conduction to obtain the thermally stationary state

solution for a surge arrester. The second treats the applicability of an envelope

equation model to simulate the heating/cooling behavior of surge arresters over

time spans which are huge with respect to the period length of the exciting electrical

signal.

While the two preceding chapters are dedicated to acquiring the model parameters

for a simulation or the simulation itself, Ch. 7 goes one step further, towards post-

8

processing. It describes the generation of a nonlinear equivalent circuit model from

the results of a numerical field simulation.

Finally, I wish to express my thanks to Siemens AG, the High-Voltage Laboratory

at TU Darmstadt and its head, professor Hinrichsen, for kindly providing some

photos and illustrations and granting me the permission to use them.

9

2 Surge Protection of Power

Transmission Lines and Field

Control

Today the protection of power transmission lines from transient overvoltages like

lightning surges or switching voltages is assured by gapless metal-oxide surge ar-

resters. Their history begins with the fortunate discovery of the nonlinear properties

of zinc oxide (ZnO) ceramics by Japanese researcher Michio Matsuoka in July 1967.

This discovery opened the way to the consequent development of ZnO varistors in

the late 1970s [72].

The appearance of the ZnO varistor had a profound effect on surge protection, as

the exploitation of their extraordinary properties made the construction of surge

arresters possible that no longer require spark gaps. Due to the technical and

economical advantages of the gapless metal-oxide surge arrester, this type of arrester

has largely supplanted its older, silicon carbide (SiC)-based predecessors, even

though some of them may still be in operation [40, 84].

The name varistor is an abbreviation of variable resistor. Varistors are also known

as voltage-dependent resistors (VDRs), which describes perfectly their characteristic

property: Their electrical conductivity depends on the applied voltage. In that

respect they act like diodes. However, they do not share the asymmetry with respect

to the direction of the applied voltage that is characteristic for diodes and, first of

all, the energy absorption capability of MOVs is much higher than that of diodes.

The combination of an extremely nonlinear electrical resistance with low leakage

currents under normal operational conditions and a relatively large heat absorption

capability are the principal factors which support their usage in surge-protective

devices (SPDs).

In its most simple form (and without paying attention to physical units), the

relationship between current I and voltage U of a varistor is described by

I = kU α , (2.1)

with material constant k and factor of nonlinearity α [21, 24]. Alternatively, α is

defined as elasticity of the function I(U):

dI dU

α= / . (2.2)

I U

10

For varistors a high (maximum) value of α is desirable and seen as a figure of

merit, because it guarantees a short transition interval between nonconducting and

conducting state. The nonlinear electric behavior will be discussed in more detail

in Sec. 3.3. The SiC varistors, which were generally used for higher voltage surge

arresters before the arrival of the MOVs, exhibit a nonlinearity coefficient in the

range of 2 to 7, whereas the values of early ZnO varistors laid in the range of 25

to 50. Clarke provides a range of 30 to 80 for commercially available varistors in

1999 [9]. Mahan et al. [69] report even an upper limit of 100 for the nonlinearity

coefficient.

The first surge arrester devices were composed of simple spark gaps. At sufficiently

high voltage differences across an air or gas gap a spark is initiated. While a spark

gap limits the occurring overvoltages, it does not ensure a fast and reliable clearing

of follow-up currents. As long as there is a voltage difference between the two

terminals of the spark gap, current continues to flow. Therefore, an alternating

current (a.c.) line-to-ground voltage signal is not interrupted until its next zero

crossing.

Following the discovery of the nonlinear properties of SiC in the 1930s, this type

of arrester has been succeeded by the valve-type arrester. The valve-type arrester

was a significant improvement over the simple spark gap designs. It combines SiC

varistors with spark gaps.

Although the elimination of the spark gaps would have been desirable to reduce

the response time of the arresters, renouncing entirely on them was not possible, as

the conductivity of the SiC varistors at normal operating voltages would have been

too high and led to too high electric losses. In consequence, the sparkover has still

been required to trigger the arrester, but the SiC varistor elements ensure a more

reliable and faster clearance of follow up currents.

Only after the development of MOVs, the gapless surge arrester became possible.

The concept of the gapless ZnO surge arrester, accompanied by experimental results,

was presented by Matsushita and General Electric in 1977/78 [59, 86].

A group of three typical (silicone-housed) station-class high-voltage MOV surge

arrresters is shown on the right-hand side of Fig. 2.1. Their design is relatively

simple. As they require no spark gaps, their only electrically active element is the

central resistor column. For practical purposes, arresters are divided into several

units. One unit of a porcelain-housed arrester is shown in Fig. 2.2. The resistor

column consists of varistors, some metallic spacers, and a supporting structure. The

spacers are used to fit the length of the column to the size of the unit, ensure an

evenly distributed contact pressure and operate as heat sinks [46]. In the figure

a porcelain-housed surge arrester with a gas-filled volume, often air, between the

central column and housing is shown. More modern arresters, especially at lower

11

Figure 2.1: Photograph with three two-unit silicone-housed surge arresters with

grading ring standing on a pedestal. Source: Siemens AG. Printed with

permission.

12

aluminum flange

metal-oxide varistor

metallic space

support structure (GRP)

air/gas volume

porcelain housing

arrester. The electrically active part is composed of the metal-oxide

varistors, metallic spacers, and its glass-reinforced plastic (GRP) support

structure. It is protected from the environment by a porcelain housing.

At top and bottom, the active part is in contact with aluminum flanges.

Between active part and the inner side of the housing a gas-filled vol-

ume is found, which is insulated from the air in the exterior. Source:

Hinrichsen [46]. Adapted with permission.

13

voltage levels, possess a silicone-rubber housing, which is often in direct contact

with the surface of the varistors by means of injection molding. Such an arrester has

no gas-filled volume between active part and interior of the housing. Therefore, no

attention has to be paid to gas pressure, humidity and leakage [46].

Basically, an arrester is supposed to behave like a voltage-controlled switch,

which is normally closed and opens temporarily during a voltage surge to allow the

transport of charges between its terminals, e.g., a transmission line and ground,

and restrict the overvoltage of the line. For most of the time, the varistors show an

insulator-like behavior and the resistive current is negligible. The measured current

between the terminals of the arresters is predominantly capacitive. This changes

during a transient voltage surge. In parallel with the voltage drop between the

terminals of the arrester, the voltage at the individual varistors increases. As soon as

the voltage reaches the strongly nonlinear region of the material characteristic, or

equivalently the U - I -characteristic, the conductivity of the MOVs increases almost

instantaneously by many orders of magnitude with response times in the nanosecond

range [43, 90]. The change from insulating to conductive state of the varistors

enables the electric charges to flow through the arrester towards ground, thereby

restricting the extent of the possible overvoltage. The overvoltage is essentially kept

from exceeding a so-called clamping voltage.

In practice, the surge arrester is designed to be operable under the maximum

root mean squared (rms) voltage of a system Us , which is higher than the nominal

system voltage. Voltage Us is defined as the highest value of voltage which occurs

under normal operating conditions at anyptime and at any point of the system [50].

The peak system-to-ground voltage (Us / 3) of such a three-phase system, plus an

additional margin of at least 5 % for eventual harmonics, has to be smaller than the

so-called maximum continuous operating voltage (MCOV) Uc of the arrester, i.e.,

the rms voltage at which the arrester can be operated continuously, demonstrated

by passing standardized test procedures defined by IEC and IEEE of 30 minutes

length [45]. Up to the MCOV of the arrester Uc , conductivity has to stay low to limit

the heat generation by electric losses. Voltages that are higher than Uc are supposed

to never occur, except for short periods during transient voltage surges. In such

an event, the arrester should be highly conductive to keep the overvoltage as low

as possible, but revert without delay to its original non-conductive state, once the

surge has passed.

Ideally, the drop of the electric potential is equally distributed among the varistors

and all of them experience the same load. However, the capacitive coupling between

the varistors and ground depends on their location in the arrester because of

stray capacitances between the resistor column and ground or nearby objects. The

installation of one or several grading rings is intended to create a more homogeneous

14

field distribution along the resistor column. Even so, at elevated voltages, e.g., the

test voltage for the accelerated aging procedure according to norm IEC 60099-4 [51],

or for larger surge arresters the distribution becomes critical to the thermal stability

of the system [39]. This problem is treated in more detail in the sections about the

calculation of an equilibrium temperature distribution (Sec. 6.1) and the possible

use of an envelope equation model to study thermal stability (Sec. 6.2).

Besides the question of thermal stability of a surge arrester after a surge, it is

important to understand, how a surge arrester reacts during fast transients like

lightning strokes or lightning test impulses. This problem has been treated by the

author [17], but is excluded from this work.

In the last few years, varistors have become available in the form of microvaristor

powder which opens up new fields for the application of ZnO varistor materials.

In particular, the possibility of using microvaristors for nonlinear resistive field or

stress control has attracted some interest.

In high-voltage engineering the stress caused by tangential electric fields at

surfaces is often a critical design parameter, for example, for cable accessories,

generator bar insulations or the fringes of a plate capacitor, as this causes discharges

that damage the respective device.

To reduce the likelihood of discharges, one tries to influence the distribution of

the electric fields by field grading. This can be done in various ways. Küchler [61]

describes five types of field grading. Changing the geometry to obtain more favorable,

rounded profiles, is known as geometric field grading. This approach is very common,

but it significantly increases the size of the respective device. A popular and effective

alternative, though not against direct current (d.c.) stress, is the so-called capacitive

field grading, whereby several sheets work as a series of capacitors to modify the

field in the desired way.

Field grading is also achieved by placing materials with appropriately chosen

properties onto the relevant surface, which alters the field distribution. Refractive

stress control is achieved by using materials with relatively high relative permittivity,

whereas resistive stress control can be realized with materials of moderately high

electric resistance, if their resistance is correctly chosen under consideration of the

capacitance between material and insulated conductor.

When the material is not simply resistive, but possesses a nonlinear, field-

dependent conductivity, one speaks of nonlinear (resistive) stress control. As the

increase of local electric field strength is slowed by the increase in conductivity,

nonlinear resistive stress control results in more homogeneous field distributions

and lower peak field strengths inside the nonlinear material. It makes more compact

designs possible than those which can be achieved with either geometric or capaci-

15

tive stress control. However, the higher conductivity leads also to higher Joulean

heat losses.

Nonlinear resistive field control with ZnO microvaristors became an option about

a decade ago, as microvaristor powders had become commercially available. Every

particle of the powder is by itself a small varistor. Apart from a smaller grain size

the structure is not fundamentally different from that of normal varistors, which

will be described in Sec. 3.3 [19, 20].

Differences in the doping and sintering process allow for the production of

microvaristor powders with largely different switching voltage and degree of nonlin-

earity.

Microvaristors can be used as functional fillers by embedding them in another

material, e.g., a polymer material like silicone rubber. The electric properties of such

a composite material are mainly determined by the type of microvaristor powder

and to a lesser extent by the volume content of the microvaristors in the composite

material. This offers an opportunity to tune the material properties to accommodate

to specific needs.

For medium-voltage cable accessories nonlinear resistive field control is already

commonly used, because of the achievable compactness of the devices and costs.

However, silicon carbide or carbon black, an industrially-produced form of soot,

are used as functional filling materials at present. The deficient reproducibility of

the composite material properties constitutes an important disadvantage. For these

materials the filler-content level has to be close to the percolation threshold, where

the material properties are very sensitive [104]. If microvaristors are used instead,

the material properties are more reproducible [30].

Beyond being a mere replacement for other materials in existing applications,

microvaristors promise the opportunity to introduce nonlinear stress control to

further areas, for which the nonlinearity of silicon carbide or carbon black has

been insufficient. For example, until now nonlinear resistive stress control of cable

accessories had been restricted to the medium voltage level, because for high-voltage

(HV) cable accessories the heat generation due to Joulean losses would have been

unacceptable. Thanks to their stronger degree of nonlinearity and high resistivity

below their switching voltage, it may be possible to build more compact accessories,

while heat development remains under control.

16

3 Physical Background

3.1 Maxwell’s Equations

almost all phenomena related to electromagnetism apart from forces and quantum

dynamical effects.

In his famous work A Treatise on Electricity and Magnetism of 1873, James Clerk

Maxwell presents a coherent mathematical system that provides an explanation for

all observed electric and magnetic phenomena and is generally known as Maxwell’s

equations. Apart from the until then unknown displacement current density, the

equations constitute no completely new discovery, but they unite the existing know-

ledge contributed by various scientists into one unified mathematical model. Most

notably, they predict the existence of electromagnetic waves propagating with a

velocity similar to that of light, a prediction that was experimentally proved correct

by H. Hertz in the late 1880s [73, 89]. Thus, light is a form of electromagnetic

wave.

Today Maxwell’s equations are usually presented as a set of four equations plus

some auxiliary equations, which describe the influence of matter on the electromag-

H R

Gauss’s law div D = % ∂Ω

D · dA = % dΩ

H Ω

Gauss’s law for magnetism div B = 0 B · dA = 0

∂B

∂ΩH R

Faraday’s law of induction curl E = − ∂t E · ds = − ∂∂ Bt · dA

∂D

H

∂A RA

Ampère’s circuital law curl H = J + ∂t H · ds = − (J+ ∂∂ Dt )· dA

∂A A

17

netic fields. In Table 3.1 the equations are shown in differential and their equivalent

integral forms (as given by Klingbeil [58]).

In these equations D, B, E, H and J are vector fields. The fields D and E are

electric field strength and electric displacement flux density, whereas H and B are

magnetic field strength and flux density. The vector J is known as current density.

the fields with matter. The vector fields in Maxwell’s equations are not independent

from each other. Their interaction depends on the environment.

In the case of vacuum, the coupling is described by the following set of constitutive

equations:

D = "0 E, (3.1)

B = µ0 H, (3.2)

J = 0. (3.3)

The electric field strength E is proportional to the dielectric flux density D. The

proportionality factor is often called electric constant and given the symbol "0 .

Similarly, the so-called permeability of free space µ0 describes the proportionality

between magnetic flux density B and field strength H.

For linear and isotropic media, one introduces material-dependent constant and

scalar relative permittivity "r , relative permeability µr and electrical conductivity σ.

The latter describes the relationship between current density J and electric field

strength E known as Ohm’s law.

B = µr µ0 H, (3.5)

J = σE. (3.6)

field interactions with smoothed out fields [26]. They give an aggregate over

volumes that are large relative to the size of individual atoms or molecules.

The relative permittivity is a consequence of the polarization of atoms and

molecules. The dipole moments of the individual molecules pn add up to a macro-

scopically averaged electric dipole moment, known as polarization density P [55]:

¬X ¶

P(r, t) = pn δ(r − rn ) , (3.7)

n

18

The polarization density P (and higher-order moments) influences the field

density D:

D = "0 E + P. (3.8)

For most materials and weaker field strengths, the polarization is proportional to

field strength E and directed into the same direction. Therefore, the influence of

polarization can be subsumed into a constant relative permittivity

Sometimes the relationship between the field quantities is not sufficiently well

described by a constant proportionality factor, e.g., in the case of magnetic field

strength H and flux density B for ferromagnetic materials. This case is of particular

practical importance, because the hysteretic behavior of ferromagnetic materials is

decisive for the functioning of most electric machines. However, this work does not

treat problems related to permeability.

The phenomenon of nonlinear conductivity σ and especially permittivity "r is of

somewhat less importance, but it is essential for the discussed varistor materials.

The field-strength dependence of the conductivity of a varistor is evident, as this is

what defines a varistor, while the nonlinearity of relative permittivity "r of a varistor

is shown to be of importance in Sec. 5.4.3.

σ = σ(E), (3.10)

"r = "r (E). (3.11)

Other cases, for which the coupling between the fields is not linear, are dispersive

and anisotropic media, the presence of impressed currents or materials with observ-

able hysteresis effects [41]. In the case of inert materials a material property, e.g.,

electric conductivity, depends not only on present field values, but also on the past

evolution of the field. This case will become the subject of Sec. 5.5.

Inert material behavior is a common phenomenon and the observation of dis-

persive, i.e., frequency-dependent, material behavior is directly connected with

the delayed reaction of a material with respect to the applied fields. In his book

about electromagnetics, Balanis [2] even presents the constitutive equations starting

directly from

Jr (t) = σ̂(t) ∗ E(t), (3.13)

19

where ∗ is used as symbol for the convolution operation. This convolution incorpo-

rates the potentially inert behavior of a medium. Present flux and current depend

on past electric fields.

Since convolution in time-domain is equivalent to multiplication in frequency-

domain, the constitutive equations in frequency domain are

Jr ( jω) = σ̂( jω)E( jω). (3.15)

Therefore, the above constitutive equations describe materials that are dispersive

or frequency-dependent due to the delayed influence of the excitation.

This delayed response is known as relaxation. For D it can be written in the form

Z ∞

D(t) = "0 "∞ E(t) + Φ̇(τ)E(t − τ) dτ , (3.16)

0

which splits the current flux density into two parts [28]. The first part is the

contribution of the momentary electrical field and the second part describes the

influence of past field values.

In the most simple case, Φ(t) is proportional to an exponential function with

relaxation time constant τ0 :

This approach is known as Debye relaxation. The exponential decay represents the

reasonable expectation of a monotonously diminishing influence of past field values.

Other models for dispersive material behavior include the Cole-Cole and Cole-

Davidson relaxation model as well as their generalization, the Havriliak-Negami

relaxation model

"s − "∞

"( jω) = "∞ + β

. (3.18)

(1 + ( jωτ0 )α )

β

Williams-Watts relaxation model (with Φ(t) = e−(t/τ0 ) ) and also the superposition

of multiple Debye relaxation models.

In Sec. 5.5 the characterization of field-strength materials with inert behavior will

be discussed introducing another generalization of the Debye model to field-strength

dependent materials.

20

3.1.3 Electro-Quasistatic Approximation

It is not always necessary to solve the full set of Maxwell’s equations. Sometimes

it is possible to neglect some part of the equations that has only a minor influence

on the solution. Whenever this is possible, the respective problem is likely much

easier to solve. We are interested in a case which we call electro-quasistatics

(EQS) and which is applicable to many problems of high-voltage engineering and

micro-electronics [83].

In electro-quasistatics, the magnetic fields are either negligibly small or they

change only slowly in time. In either case, the term ∂∂ Bt in Faraday’s law becomes

negligible. This implies that there are no eddy currents and that the electric field

strength E can be regarded as irrotational. Given that E is irrotational, it can be

represented as the gradient of a scalar potential ϕ :

E = − grad ϕ. (3.19)

∂D

div curl H = div J + (3.20)

∂t

∂

= div σE + ("r "0 E) (3.21)

∂t

∂

= − div (σ grad ϕ + ("r "0 grad ϕ) . (3.22)

∂t

∂

div σ grad ϕ + ("r "0 grad ϕ) = 0, (3.23)

∂t

or:

∂

div (σ grad ϕ) + div ("r "0 grad ϕ) = 0. (3.24)

∂t

It is this last equation that is being solved for EQS problems.

Applicability of quasistatics

While the equation used for electroquasistatics simulations has been derived in the

preceding section, it is not obvious, in which cases it can be applied. One needs

some indicators to decide, if it is appropriate for a specific problem or not.

21

At first, one can look at the energy of the electric and magnetic fields. If the

energy stored in the magnetic field wm is much smaller than the energy stored in

the electric field we , i.e.,

wm we , (3.25)

then the choice of electro-quasistatics is probably more appropriate than magneto-

quasistatics [100]. The energies are given by [44]:

ZB

wm = H · dB, (3.26)

0

Z D

we = E · dD. (3.27)

0

If the inverse holds true and the magnetic energy is much larger, then the problem

tends to be magneto-quasistatic.

Time constants are another, more relevant help to decide whether a problem is

magneto- or electro-quasistatic [91]. There are three time constants that have to be

considered:

• characteristic time τ,

The first time constant, the characteristic time τ determines whether the full

Maxwell system has to be solved or if a simplification is possible. The other two

time constants decide over the kind of simplification that can be made.

If the characteristic time τ does not fulfill the condition

L

τ , (3.28)

c

then neither the electro-quasistatic nor the magneto-quasistatic simplification is

valid. The characteristic time τ of the equation above is defined as

1 1

τ= = , (3.29)

ω 2π f

where f is the frequency of a sine-like excitation. L is the typical length of the

considered geometry and c is the wave propagation speed in the medium. Essentially,

equation 3.28 requires that phase differences in the model are negligible, since:

L

τ =⇒ λ L, (3.30)

c

22

with wave length λ = c/ f .

The presented inequality relations follow from an estimation of the maximum

relative error for the electric field strength E, when the term ∂∂t B is neglected in

electro-quasistatics, and of the error for the magnetic field strength H due to the

disregard of ∂∂t D in magneto-quasistatics. U. van Rienen [99] deduces the condition

by substituting the time derivatives by a factor 1/τ and the space derivatives by 1/L

in Faraday’s law of induction and Ampère’s circuital law.

The maximum estimated induced electric field strength E due to magnetic field

strength H is

µr µ0 L

E= H. (3.31)

τ

"r "0 L

H= E. (3.32)

τ

simplification is made, can be estimated by inserting equation 3.32 into equation

3.31. The retroactive change of field strength Er becomes

µr µ0 L "r "0 L L2

Er = E = 2 2 E. (3.33)

τ τ c τ

3.28 arises as a direct consequence. For magneto-quasistatics, the same argument is

made, except that equation 3.31 is inserted into equation 3.32 and Hr H has to

be fulfilled.

The characteristic time constant τ helps to decide, if a quasistatic approach is valid.

But it gives no answer to the question, if a problem is electro- or magneto-quasistatic.

There are two other time constants that have to be considered. The first one is

the dielectric relaxation time τrelax [65],

"

τrelax = , (3.34)

σ

23

The dielectric relaxation time is a consequence of the continuity equation, which

is itself derived from Ampère’s and Gauss’s law. In a homogeneous medium, one

has:

∂

0 = div J + D (3.36)

∂t

∂

= div J + div D (3.37)

∂t

σ

∂

= div D + div D (3.38)

" ∂t

σ ∂

= div D + div D (3.39)

" ∂t

σ ∂%

= %+ . (3.40)

" ∂t

By solving Eq. 3.40, one obtains Eq. 3.35 with dielectric relaxation time τrelax .

The other time constant is called magnetic diffusion time τdiff ,

needs to diffuse into a previously field-free domain of characteristic length L , e.g., a

cube with side length L .

Starting from Ampère’s circuital law with displacement currents and assuming a

homogeneous medium, the magnetic diffusion time can be derived as follows:

curl curl H = curl (σE) , (3.43)

1

curl curl B = σ curl E, (3.44)

µ

1 ∂B

(grad (div B) − ∆ B) = −σ , (3.45)

µ ∂t

∂B

∆ B = µσ . (3.46)

∂t

Faraday’s law of induction and Gauss’ law for magnetism are invoked in equations

3.44 and 3.45.

As previously done for characteristic time τ, the space and time derivative opera-

tors are substituted by factors 1/L and 1/τdiff . It follows directly that τdiff = µσL 2 .

24

Finally, one may observe that the characteristic time τ, the dielectric relaxation

time τrelax and the magnetic diffusion time τdiff are not mutually independent. The

characteristic time τ is the geometric average of the other two time constants:

p

τ= τrelax τdiff . (3.47)

Therefore, the dielectric relaxation time has to be much larger than the magnetic

diffusion time for an electro-quasistatic problem:

from a body without work and which is caused by temperature differences [82].

The heat conduction equation, which describes the transport of heat in solids,

was formulated and presented to the members of the Institut National at Paris by

Jean Baptiste Joseph Fourier in 1807 [32]. Thus, it precedes Maxwell’s treatise,

which laid the foundations of electromagnetic theory, by 60 years.

However, it took more than fifteen years until Fourier’s theory became widely ac-

cepted. In 1808 Poisson presented Fourier’s results and the relevant heat conduction

equation in an article in the bulletin of the Societé philomathique, but the theory

was rejected by its reviewers, among them Lagrange and Laplace, due to ostensible

gaps in the application of trigonometric series. Only after the publication of Fourier’s

book Théorie analytique de la chaleur in 1822 his theory gained momentum [31, 78].

In vector calculus notation, the equation can be presented as 1

∂ϑ

cp % − div (λ grad ϑ) = 0. (3.49)

∂t

capacity per mass of a medium, % is the mass density, and λ is known as thermal

conductivity. The product of c p and % may be combined into volumetric heat

1

In general, the

change of heat capacity is neglected. although, in principle, the first term should

be ∂∂ t c p %ϑ .

25

capacity cv = c p % . Introducing a diffusivity α = cλv the equation can be simplified

even further [37].

Physically, the equation represents an energy balance. The term c p % ∂∂ ϑt describes

the rate of change of the energy stored per volume, while the second term stands

for the heat flowing into or out of the volume. Fourier’s law posits that the heat flux

density q̇ is proportional to the gradient of temperature:

q̇ = −λ grad ϑ. (3.50)

Therefore, div q̇ is the net heat flow out of an infinitesimally small volume. An

eventual increase or decrease of the local thermal energy has to be accompanied

by the net inflow or outflow of energy through the boundary of the volume, as

the original version of the heat conduction equation does not consider local heat

sources.

By including local heat sources, e.g., Joulean losses, nuclear fission or exothermic

chemical reaction, the heat conduction equation becomes [12]

∂ϑ

cp % − div (λ grad ϑ) = ẇgen . (3.51)

∂t

In this equation, ẇgen designates heat generation density.

The heat conduction equation describes the propagation of heat inside of solid

materials. In fluids and gases the principal mechanism for the transportation of heat

is not conduction, but either convection or radiation. Therefore, the heat conduction

equation may not be applied in these cases. Its field of application is restricted to

solid bodies. By contrast, convection and radiation can be neglected in the context

of heat propagation in solids.

Nevertheless, convection and radiation are important boundary conditions to

heat propagation problems. Additionally, one may remark that the individual heat

transportation mechanisms are not mutually exclusive. Heat conduction, convection

and radiation can take place at the same time, even if typically one mechanism is

dominant.

Convection

Heat conduction is related to the interaction of particles in the molecular range, for

example, lattice vibrations and diffusion of electrons in solids. Macroscopically there

26

is no displacement of matter [63]. By contrast, in fluids and gases the molecules,

carrying energy in the form of enthalpy, move over much greater distances, leading

to convective heat transfer.

If this heat transfer occurs naturally, one speaks of free or natural convection.

Otherwise, if the heat transfer is sustained by technical means or other external

factors, one speaks of forced convection [53].

In principle, the determination of the velocity field of the particles and of the con-

vective heat transfer requires the solution of the Navier-Stokes equation. However,

Newton’s law of cooling is much simpler and often sufficient to describe the heat

transfer between the surface of a solid and a fluid or gas:

q̇ = h (ϑ − ϑ∞ ) . (3.52)

In this empirically derived equation ϑ is the temperature at the surface of the solid

body and ϑ∞ is the far-field temperature in the gas or fluid, i.e., the temperature

at a sufficiently large distance from the thermal boundary layer [57]. Unlike heat

conductivity λ, the heat transfer coefficient h is not exclusively material-dependent,

but it depends also on geometry, solid-surface characteristics (e.g., roughness) and

other flow structure aspects.

For special cases, values for the heat transfer coefficient h can be found in tables

with experimental results, otherwise one has to solve some boundary layer equation

or even the Navier-Stokes equations.

Radiation

Instead it increases with the fourth power, so that it tends to become the predominant

mechanism of heat transfer at higher temperatures [76, 82].

It is related to changes of molecular energy levels. When an electron drops from

a higher level to a lower one, a photon, or equivalently electromagnetic radiation, is

emitted, mostly at wavelengths between 100 nm and 1 mm. Because of its nature as

electromagnetic wave, radiative heat transfer can take place over large distances

and through vacuum, in contrast to conduction and convection.

Radiative heat transfer is not easy to model, because, as stated by Modest [76],

radiative properties are usually difficult to measure and often display erratic behavior.

For liquids and solids the properties normally depend only on a very thin surface layer,

which may vary strongly with surface preparation and often even from day to day. All

radiative properties (in particular for gases) may vary strongly with wavelength, [. . . ].

The heat flux density emitted from a gray surface is given by

q̇ = εσSB ϑ4 , (3.53)

27

with emissivity ε, Stefan-Boltzmann constant σSB and surface temperature ϑ. A

black emitter is an idealized emitter with an emissivity of one. It absorbs all radiation

and emits heat radiation of a specific frequency spectrum. Real materials have a

emissivity lower than one and are frequently modeled as so-called gray emitters. A

gray emitter is an idealized material with the spectrum of a black emitter, but which

absorbs and emits proportionally less radiation (0 < ε < 1) [1].

In practice, one is interested in the net heat flux between a body and its envi-

ronment. The net heat flux between two non-transmissive surfaces that can be

considered as grey and diffuse emitters has the following form

Q̇ = Σ12 T14 − T24 , (3.54)

with coefficient Σ12 . Non-transmissive means that the incident radiation is com-

pletely absorbed or reflected at the surface. Diffuse means that the photons are

emitted in all directions with the same probability.

The coefficient Σ12 is a function of the two emissivity values and geometry. The

value for two large parallel plates is

σSB

Σ12 = 1 1

A. (3.55)

ε1 + ε2 −1

In the case of two concentric cylinders, which will arise in Sec. 6.1, the situation

becomes still more complicated than for two parallel plates, as one has to take

into account the kind of reflection on the surfaces. One distinguishes between

specular and diffuse reflection. Specular reflection means that the incident radiation

is reflected according to the law of reflection, like light at a mirror, while ideal

diffuse reflection means, as explained above, that the radiation is reflected into all

directions equally. For two concentric cylindrical surfaces, the value for Σ12 depends

entirely on the kind of reflection of the outer surface, which determines the outcome

of the re-reflections [48]:

A1 σSB

1 1 if A2 specular

ε1 + ε2 −1

Σ12 = A1 σSB (3.56)

1 A1 1 if A2 diffuse,

ε1 + A2 ( ε2 −1)

where the index 1 is used for the inner and 2 for the outer surface. It is noteworthy

that the area A2 has no influence on Σ12 in the case of a specular surface A2 . In

practice, the real heat exchange is most often a combination of the models as a

surface is rarely perfectly specular or diffuse.

28

3.3 Metal-Oxide Varistors

extraordinarily nonlinear electrical conductivity. This behavior has been observed

for various metal oxides, for example, aluminum and titanium oxides, but the ZnO

varistor is the oldest and by far most common type of varistor.

ZnO varistors are heterogeneous ceramic materials. They consist of about 90 % or

more of ZnO. Their manufacturing process is described extensively by Einzinger [24]

and more recently by Clarke [9] or Haddad [40]. First, very fine powders of almost

pure ZnO and of doping materials like Bi2 O3 are mixed and ball milled over many

hours to obtain a homogeneous mixture. After addition of a binder material the

aqueous slurry is calcinated at about 700 ◦C. During this process the dopants dissolve

homogeneously in the solution.

Then the resulting granulate is compacted, molded into varistor shape and

sintered in a furnace at temperatures above 1000 ◦C. During the sintering process

the particles of the ZnO granulate grow together and form larger grains.

After sintering the material cools slowly and in controlled way to influence the

chemical processes that take place at the shared surface of two grains and at the

multi-grain junctions. At the boundary between two grains a thin, dopant-enriched

intergranular layer forms. This inter-grain boundary is decisive for the electrical

properties of the future varistor.

The final steps of the varistor manufacturing process include the coating to reduce

the risk of surface sparkovers and the deposition of metals, e.g., aluminum, at the

contact surfaces.

Since the beginning of research on varistors, its nonlinear electrical behavior

has been attributed to the grain boundaries [71]. While the grains, which have a

diameter in the order of 10 µm, fill a great part of the volume, they are composed

of almost pure, crystallized ZnO, which is known to be a moderately good, linear

electrical conductor (σ ≈ 100 S/m). Therefore, the explanation for the nonlinear

varistor behavior has be found at the boundaries. Notwithstanding the achieved

progress, the actual mechanism remains not completely understood [56]. Over the

years, many different mechanisms have been proposed, spanning from diffusion

theory, avalanche theory and tunneling to more refined models like double Schottky

barriers with thin layers or depletion layers with recombination of holes [40].

Measurements show that the grain boundaries are insulating up to a switch-

ing voltage of about 3.5 V. Published values for different dopant additives and

manufacturing processes lie in the range between 2 V and 4 V [38]. Furthermore,

29

the conductivity may vary between different boundaries of the same varistor, as

confirmed by Tao et al., who divided the grain junctions into "good" and "bad" [96].

Above switching voltage the boundary becomes conductive and the bulk conductivity

approaches the value of the ZnO grains.

The varistor microstructure with its grains and boundaries of different size and

electrical properties affects also the observed failure mechanisms of varistors, notably

the so-called puncturing. In that case, privileged current paths develop between

the contact surfaces of the varistor, where current has to pass through less, or less

resistive, boundaries to reach the opposite side. Visually, this failure mode can be

identified by looking for small spots on the electrodes, where the highly localized

current reaches the surface. This failure mode is associated with long-duration,

low-amplitude impulses and the melting of the Bi2 O3 additive at temperatures of

about 820 K [22, 85]. It has to be noted that this temperature may be reached at

the grain boundaries during an impulse, while the temperature of the grains is still

much lower.

If the impulse is shorter and of larger amplitude, the predominantly observed

failure mode becomes cracking. In that case, thermal stress causes the varistor to

crack, possibly bursting apart in pieces and requiring measures to protect nearby

objects or persons. A third observed reason for failure are surface flashovers. The

above mentioned coating of the surface serves to reduce the probability of this kind

of failure.

100

σ(E, ϑ) = 200 ′

1 a1

2A1 +A2 exp ϑ

a a

1−tanha2 + ϑ3 − E4

Conductivity / S/m

10−2

10−4

493 K

10−6 443 K

393 K

10−8 343 K

10−10 293 K

1 1.5 2 2.5 3

Field Strength / kV/cm

Bartkowiak [5]. The electrical conductivity varies by a factor of up to

1 · 1012 at normal ambient temperatures (≈ 293 K).

30

a) b)

c) d)

f)

C

e) g)

Figure 3.2: Some varistor equivalent circuit models: a) Harnden et al. (1972, [43]), b)

Fernández and Díaz (2001, [29]), c) Pinceti and Giannettoni (1999, [81]),

d) IEEE frequency-dependent model (1992, [52]), e) Eda two-grain model

(1989, [21]), f) Einzinger (1978, [23]), g) parallel nonlinear resistor and

capacitor.

31

The fourth and possibly most important destruction mechanism is the so-called

thermal runaway. Due to the temperature dependence of the electrical conductivity

of the varistors, after one or multiple voltage surges a surge arrester may reach a

state, from which it does not return to normal operational conditions. Instead, the

arrester continues to heat up leading to the destruction of the varistors.

microstructure. In Fig. 3.1 the electrical conductivity is shown as a function of

electric field strength and temperature. The diagram is derived from an empirical

I-U-characteristic presented by Bartkowiak [4]. While the characteristics of varistors

vary significantly, the diagram illustrates clearly the extreme nonlinearity typical for

varistors. At high field strengths the conductivity is up to 1 · 1012 times higher than

at low field strengths.

Another important property, which can be seen in Fig. 3.1, is the temperature de-

pendence. The conductivity increases, as temperature rises. Equivalently, resistivity

decreases. Materials with decreasing resistivity as temperature increases possess a

negative temperature coefficient (NTC) and are known as NTC resistors. Often, as

in the case of surge arresters, this property is undesired, as it may facilitate heat

and stability problems, known as thermal runaway. The thermal stability of surge

arresters is an important problem, which will be discussed in more detail in Sec. 6.2.

ZnO varistors possess also a surprisingly high relative permittivity. The macro-

scopic permittivity is much higher than the permittivity of its composing materials.

A varistor is largely filled with ZnO which has a relative permittivity in the range

of 8 to 10. The relative permittivity of the dopant materials does not exceed 25.

Nevertheless, measurements show that the relative permittivity lies in the range of

several hundreds. Haddad [40] or Levinson and Philipp [67] even give values in

the range of 1000 to 1600 for the calculated relative permittivity from capacitance

measurements. These surprisingly high values are only partially explained by the

structure of the material, where the charges are separated by very thin grain bound-

aries. In addition, depletion layers adjacent to the intergranular layer and trapped

electrons are cited as explanation.

Numerous equivalent circuit models have been proposed for modeling the electri-

cal behavior of ZnO varistors. Fig. 3.2 presents a non-complete list of equivalent

circuit models. The circuit elements labeled with A0 or A1 are idealized, purely

resistive voltage-dependent resistors or varistors. Model f ) of Fig. 3.2 deserves some

special attention, as it does not only include a nonlinear resistor element but also a

nonlinear capacitor. The results of the (micro-)varistor characterization in Sec. 5.4.3

32

Volumetric Heat Capacity / J/(cm3 ◦C)

3.5 3.5 · 106

cv = c0 + α c T

3.0 a 3 · 106

e b

2.0 2 · 106

c d

1.5 1.5 · 106

Temperature / ◦C

700 7 · 108

αc

E/V = c0 T + 2 T2

600 6 · 108

Energy Density / J/cm3

500 5 · 108

b 4 · 108

400

300 a d 3 · 108

e

200 2 · 108

c

100 1 · 108

0 0

0 50 100 150 200

Temperature / ◦C

Figure 3.3: The two diagrams show the nonlinear volumetric heat capacity of some

surge arrester materials as given by Lat [64] and the required energy to

reach another temperature starting from 0 ◦C. (a: ZnO; b: porcelain, c:

epoxy with SiO2 ; d/e: grey/black EPR rubber.)

33

will show that permittivity is indeed field-strength dependent. Yet, an explicitly

nonlinear capacitor is not necessary for a macroscopically nonlinear, field-dependent

permittivity. At the terminals the imaginary part of the impedance is nonlinear for

the other models, too. The finite element model according to the description in

Sec. 4.2 corresponds to a network of parallel, nonlinear resistors and conductors.

The coupling between the surfaces of the individual elements can be represented

by a resistor and capacitor in parallel as shown in the model of Fig. 3.2 g). This

model is necessarily macroscopic, as the geometric size of the elements is many

times larger than the ZnO grains.

While the thermal properties of ZnO varistors are not as extraordinary as the elec-

trical conductivity, they are nonetheless noteworthy and relevant for the functioning

of varistors and surge arresters. ZnO possesses a very high thermal conductivity.

It is for this reason that ZnO is used as additive in rubber tires [56]. A study on

the thermal conductivity of ZnO varistors is found in Barrado et al. [3]. Another

property is the high heat absorption capability of ZnO varistors. Both properties are

beneficial for the operation of MOV surge arresters. In Fig. 3.3 the heat capacity

of a ZnO varistor according to Lat [64] is shown along to those of various other

surge arrester materials. It can be seen that ZnO exhibits a comparatively high heat

capacity and absorbs, in consequence, significantly more thermal energy for the

same effective increase in temperature.

34

4 Finite Element Method

4.1 Introduction to the Finite Element Method

The finite element method is a popular approach for the numerical solution of

partial differential equations. Other well-known methods are the finite differences

method [94, 103], finite integration technique [101], finite volume method [27] or

boundary element method [34].

The name finite element method was introduced by Clough in 1960 [10]. How-

ever, the method can be traced back to the 1940s and the pioneering work of

Hrennikoff, McHenry, Courant and others, who began to apply the method for

structural analysis problems of airplanes [11, 49, 74].

Notwithstanding its successful application and the advances made by Argyris,

Turner and many others, the method remained for many years a tool exclusively used

in mechanical engineering. In 1965 its wider scope was remarked by Zinkiewicz

[108], who wrote that the method could be interpreted in terms of variational

procedures and that the minimization of the total potential energy functional leads

naturally to its extension to other boundary value field problems.

He pointed explicitly to problems of heat conduction, for which Wilson and

Nickell [102] applied the method in the following year. The first known application

of the finite elements method in electromagnetics dates to 1969 and the analysis of

a waveguide by Silvester [88].1

Today the finite element method is applied in many areas of physics and beyond.

Its application requires the execution of several steps:

which a finite element formulation can be applied.

locally defined basis and weighting functions.

3. The entries of the local matrix and force vector are calculated for each

element. A global system matrix and right-hand side vector, often called force

vector, are assembled from the local data.

1

Overviews over important contributions to the finite element method can be found in an article

by Zinkiewicz [107] and several books, e.g., [13, 106].

35

4. After the assembly of a global system matrix and force vector, the system is

solved.

In many cases, some of the previous steps are executed repeatedly, e.g., for a

transient analysis or for mesh refinement.

When applying the finite element method, one is interested in some unknown

function u(x), which is the solution of a differential equation. Instead of the

analytical solution, one searches a vector ũ(x) to approximate the true solution:

X

u(x) ≈ ũ(x) = ũi Ni (x), (4.1)

i

Depending on the choice of approach for obtaining ũ(x), one distinguishes be-

tween variational and weighted residual method. The variational method is based

on the minimization or maximization of a scalar functional Π, which reaches its ex-

tremal value for the solution of the considered differential equation. The functional

Π is described by an integral of the form

Z Z

∂u ∂u

Π= F u, , . . . dΩ + E u, , . . . dΓ . (4.2)

Ω

∂ xi Γ

∂ xi

with differential operators F (. . .) and E(. . .). Often the functional has a physical

meaning, representing energy or another relevant property [106].

After insertion of Eq. 4.1 into the functional Π, the functional is differentiated

with respect to all components of vector ũ and the result is set to zero to find a

stationary point:

∂Π

∂ ũ1

∂Π ..

= . = 0. (4.3)

∂ ũ ∂Π

∂ ũn

If neither the function nor its derivatives occur in powers exceeding two, ∂∂ Πũ = 0

provides a linear system of equations that can be solved to get a finite elements

solution of the differential equation.

Alternatively, a weighted residual method can be applied. Since this requires

neither knowledge nor existence of a variational formulation, it is much more

generally applicable. If the differential equation can be written as

Lu = f , (4.4)

36

then the residual of its numerical solution ũ is

r = L ũ − f . (4.5)

For the exact solution, the residual is zero everywhere. Thus, integration over an

arbitrary domain Ω gives

Z

r dΩ = 0. (4.6)

Ω

Obviously, the result is not affected, when the residual is weighted by a multiplicative

factor w:

Z

r w dΩ = 0. (4.7)

Ω

If the exact solution u is replaced by the numerical solution ũ, then the residual

is no longer zero, except for some points in space or in exceptional cases. Still, one

can require that the residual becomes zero in the weak sense of equation 4.7 for a

set of specified functions w i . In that case, the numerical solution is no longer exact,

but one may expect that it is still similar to the true solution. One speaks of a weak

form.

Depending on the choice of functions (or distributions) w i , several well-known

methods can be classified as weighted residual methods [6, 106]:

for ũ and weighting functions (w i = Ni ).

weighting function, i.e., the residual at certain points must be zero.

mains. In the j -th subdomain w j assumes value one and is zero elsewhere.

popular choice and is used for the simulations in the later sections. A great advantage

of this method is that the mass and stiffness matrices are symmetric, which is an

important property for iterative and direct solvers.

The finite element standard procedure is based on the individual treatment of

each element. Assuming that a weighted residual method has been chosen, so that

37

one can start from Eq. 4.7 and that the residuals are linear in ũi , the condition

related to the j -th weighting function can be written as

Z

0= r(ũ)w j dΩ (4.8)

Ω

XZ

= ri (ũi )w j dΩ (4.9)

i Ω

XXZ

= ri (ũi )w j dΩk , (4.10)

k i Ωk

with Ωk denoting the subdomain of the k-th element. Therefore, the integrals have

to be evaluated for all combinations of (i, j, k). However, if the basis and weighting

functions are zero on most elements, the computational effort is much lower. For

this reason, one chooses functions that are nonzero on only a few neighboring

elements. Thanks to their compactness, only a few combinations of i and j have to

be evaluated for any element k. Generally, the functions ũi and w j are defined as

polynomials inside the elements, in which they are nonzero.

Typically, the evaluation of the integrals in Eq. 4.10 is executed locally. This

means that the global system of equations is not assembled directly, but that one

iterates through all elements. For every element, a local element matrix Ak and

vector b k is constructed as first step. This local matrix and vector depend only on

the local basis and weighting functions and are of much lower dimension than the

global system. After construction of Ak and b k , the entries are inserted into the

global system matrix A and force vector b requiring a translation of the local row

and column number to the global ones. The resulting global system will be sparse

due to the limited coupling between the basis and weighting functions.

The following explanations for the cases of electro-quasistatics and heat con-

duction are presented in a way that avoids the discussion of all aspects related to

individual elements. However, one has to be aware that the mass and stiffness

matrices, as well as other matrices and vectors, are at first built locally, i.e., for

individual elements, but that the equations are solved at the global level.

∂

div (σ grad ϕ) + div ("r "0 grad ϕ) = 0. (4.11)

∂t

38

Therefore, a Galerkin finite element solution is obtained by solving

Z Z

∂

div (σ grad ϕ) w i dΩ + div ("r "0 grad ϕ) w i dΩ = 0, (4.12)

Ω Ω

∂t

theorem, the first integral, which is related to resistive currents, becomes

Z

div (σ grad ϕ)w i dΩ

Ω

Z Z

= div (σ grad ϕw i ) dΩ − σ (grad ϕ) T grad w i dΩ (4.13)

Ω Ω

Z Z

∂ϕ

= σ w i dΓ − σ (grad ϕ) T grad w i dΩ. (4.14)

Γ

∂n Ω

Γ denotes the surface of space Ω and the corresponding integral represents the

weighted current flowing through the surface. This integral is of importance for the

imposition of boundary conditions, but for the moment, it may be ignored. It will

be discussed in Sec. 4.2.3. By choosing the Bubnov-Galerkin

P approach (w i = Ni )

and substituting ϕ by the numerical potential ϕ̃ = j ϕ̃ j N j , the second integral on

the right-hand side becomes

X Z

T

ϕ̃ j σ grad N j grad Ni dΩ. (4.15)

j Ω

As a result of the preceding operations, one spatial differentiation has been trans-

ferred from the basis function N j to weighting function Ni , so that the requirements

on the continuity of the basis functions are reduced. The evaluation of the integral

for a combination of i and j gives the coefficient of the so-called stiffness matrix Kσ

in row i and column j . The name stiffness matrix and the letter K are a testimony

to the origins of the finite element method in structural dynamics. Because of the

interchangeability of Ni and N j , the matrix is symmetric.

The second integral of Eq. 4.12, which is related to displacement currents, is

transformed in the same way as the first integral. However, we will assume that

the relative permittivity "r is a function of the absolute value of field strength, i.e.,

"r = "r (|E|). This is done with respect to the observations that will be discussed in

Ch. 5 and leads to

39

Z

∂

div ("r "0 grad ϕ) w i dΩ (4.16)

Ω

∂t

Z Z T

∂ ∂

= div ("r "0 grad ϕ) w i dΩ − ("r "0 grad ϕ) grad w i dΩ

∂t ∂t

ZΩ ZΩ T

∂ ∂

= n· ("r "0 grad ϕ) w i dΓ − "0 ("r (|E (t)|) grad ϕ) grad w i dΩ.

Γ

∂t Ω

∂t

Again, the surface integral may be ignored for the moment. The stiffness matrix

K" is obtained from the remaining term in the same way that Kσ was obtained.

Neglecting the boundary integrals, the electro-quasistatics equation can be written

in the form

∂ ϕ̃

K" + Kσ ϕ̃ = 0, (4.17)

∂t

or alternatively

∂ ϕ̃

= −Kσ −1 Kσ ϕ̃. (4.18)

∂t

Both electric conductivity σ and relative permittivity "r are allowed to be nonlinear.

This assumption is made to correctly simulate the behavior of (micro-)varistor

materials. It has important implications on the approach for time integration.

First, the problems involving the simulation of varistor behavior are highly non-

linear. The conductivity of a varistor increases by many orders of magnitude, once

a threshold value, usually called switching voltage, is reached. Relatively small

changes of field strength result in a considerable change of the conductivity value.

Therefore, the problem becomes stiff and one has to use an implicit time-integration

scheme.

The meaning of stiffness is not precisely defined and there exists no universally

accepted definition. Deuflhard [18] clarifies:

"There are initial value problems, for which [. . . ] explicit one-step meth-

ods require a too small [length of time step] τ∆ and thus too high an

effort [. . . ]. Such initial value problems are called stiff in the literature,

all other initial value problems are called nonstiff. This definition has to

40

remain necessarily vague, because an assessment like too small can only

be done pragmatically by estimating the computational costs and then

taking the decision for or against a class of methods. Here, it is not the

problem class that leads to the choice of a method, but a class of methods

classifies the problems!" 2

Eq. 4.18 has this form. However, it requires knowledge of matrix K"−1 . This matrix

is a priori unknown and has to be obtained by inversion of K" . As the finite element

stiffness matrix K" is nondiagonal, one incurs high computational costs by this

operation. Furthermore, the sparsity pattern is not maintained. Unless the matrix

can be reused frequently, one prefers to avoid the matrix inversion. In our case

matrix K" is needed many times during time integration, at every time step, but it is

not constant. K" depends on field strength and has to be recalculated frequently, as

the electric fields change in time.

∂ ϕ̃

To avoid this problem, one can opt for solving first for K" ∂ t and then solve the

n+1

resulting linear system of equations to get ϕ̃ :

K" n n

y == ϕ̃ + f (Kσ n , ϕ̃ n , Kσ n+1 , ϕ̃ n+1 ) (4.20)

∆t

K" n+1 n+1

ϕ̃ =y (4.21)

∆t

The price of avoiding the matrix inversion is the need for solving an additional

system of linear equations.

There exist several methods for the implicit integration in time, which can be

divided into one- and multi-step methods. Runge-Kutta methods are a prominent

example of one-step methods, while Rosenbrock-Wanner methods are a popular

choice among multi-step methods.

In our simulations the θ -Method was used instead, as singly-diagonally implicit

Runge-Kutta (method) (SDIRK) methods, which belong to the family of Runge-Kutta

2

Translation by the author. Original text: "Es gibt nun Anfangswertprobleme, für welche die [. . . ]

expliziten Einschrittverfahren zu kleine τ∆ und damit zu großen Aufwand benötigen [. . . ]. Solche

Anfangswertprobleme werden in der Literatur steif genannt, alle anderen Anfangswertprobleme

heißen nichtsteif. Diese Definition muss notwendigerweise vage bleiben, da Einschätzungen wie

zu klein nur pragmatisch erfolgen können, indem der Rechenaufwand bewertet wird und sodann

die Entscheidung für oder wider die Verfahrensklasse getroffen wird. Hier führt also nicht eine

Problemklasse zur Wahl von Verfahren, sondern eine Verfahrensklasse klassifiziert die Probleme!"

41

methods, with or without Jacobian matrix delivered at best comparable simulation

times. The θ -method can be written as

ϕ̃ n+1 − ϕ̃ n

= θ F (t n+1 , ϕ̃ n+1 ) + (1 − θ )F (t n , ϕ̃ n ), (4.22)

∆t

or

In this equation ϕ̃ n is the present solution and ϕ̃ n+1 the subsequent solution after

advancing in time by a time step of length ∆t . A comparison with equation 4.18

shows that in our case

Several choices of θ lead to methods that are known under special names. In

particular, θ = 0 is the explicit or forward Euler method and θ = 1 is known as

implicit or backward Euler method. The second-order accurate Crank-Nicolson

method is obtained by the choice θ = 12 and the value θ = 23 is related to the name

Galerkin [8].

A slightly different equation is obtained from Eq. 4.12 by integration over one

time step:

Z t n+1 Z Z t n+1

div (σ grad ϕ̃) Ni dΩ dt + div ("r "0 grad ϕ̃) Ni dΩ = 0. (4.25)

tn Ω Ω tn

This leads to

Z t n+1

K"n+1 ϕ̃ n+1 − K"n ϕ̃ n =− Kσ (t)ϕ̃(t)dt (4.26)

tn

≈ −∆t θ Kσn+1 ϕ̃ n+1

+ (1 − θ ) Kσn ϕ̃ n

. (4.27)

K"n+1 + θ ∆t Kσn+1 ϕ̃ n+1 = K"n − (1 − θ )∆t Kσn ϕ̃ n . (4.28)

42

4.2.3 Boundary Conditions

ing their importance. The two basic boundary conditions for electro-quasistatics

are the Dirichlet and Neumann boundary conditions.. By Senior and Volakis [87],

these conditions are labeled the classical boundary conditions for scalar fields. As

the Robin boundary condition is not very different and is needed anyway in the

context of heat conduction, it will be included in the following presentation of

electric boundary conditions.

depends linearly on potential ϕ̃ .

These three types of boundary conditions are sufficient for all the problems,

which will be discussed in the following. However, it has to be stated that this

is no exhaustive list of boundary conditions in computational electromagnetics.

For example, corresponding boundary conditions are applied with respect to the

magnetic field. In addition, several boundary conditions have been developed to

approximate the field close to the surface of materials, which remain intentionally

outside of the computational domain, so-called impedance boundary conditions

(IBCs). Another extremely important type of boundary condition are non-reflective

or absorbing boundary conditions (ABCs), including the perfectly matched layer

(PML). They are used to terminate the computational domain for problems of waves

that propagate in unbounded regions [95].

The application of Dirichlet boundary conditions is straightforward. The potential

ϕ̃ is known in some parts of the domain and the corresponding degrees of freedom

can be eliminated from the system of equations. Sometimes one prefers not to

change the dimension of the system and keeps the entry on the diagonal of the

system matrix and sets the value on the right-hand side to the value needed to

impose the Dirichlet condition.

The Neumann and Robin boundary conditions enter the mathematical model

through the boundary integrals that result from the application of the divergence

theorem in Sec. 4.2.1. Neglecting the surface integral in Eq. 4.14 is equivalent to

imposing the so-called homogeneous Neumann boundary condition, whereby no

current flows through the boundary.

In principle, a Neumann boundary condition prescribes the normal derivative of

the unknown field on the boundary, i.e., the field strength in normal direction, but

43

this is almost equivalent to the practically more important imposition of a current

density. Both conditions can be implemented through evaluation of the boundary

integrals:

Z Z

∂ϕ

σ w i dΓ = − σ (ET · n) w i dΓ (4.29)

Γ

∂n

ZΓ

=− (JTr · n) w i dΓ . (4.30)

Γ

Here, E designates the electric field strength and Jr the resistive current density.

The products ET · n and JTr · n are the normal components of the electric field and

current density with respect to the boundary.

If the displacement current boundary integrals are disregarded, Jr is the total

current flowing through the boundary. However, formally a similar expression arises

for the displacement current density Jc :

Z Z

∂ ∂ϕ ∂

"r "0 w i dΓ = − "r "0 ET · n w i dΓ (4.31)

Γ

∂t ∂n Γ

∂t

Z

=− (Jc · n) w i dΓ . (4.32)

Γ

This term will be disregarded and we will assume that a current is only imposed,

where its capacitive part is negligible, as is the case for the examined models. By

evaluation of the integrals for all i with weighting functions that are non-zero on the

concerned boundary, one obtains a vector b, which will modify only the right-hand

side of the global linear system.

The Robin boundary condition is similar to the Neumann boundary condition,

but it will also influence the left-hand side as it depends on the local field. A Robin

boundary condition can be written in the form:

∂u

= αu + β. (4.33)

∂n

This is inserted into the boundary integral:

Z Z

∂ ϕ̃

σ Ni dΓ = σ (αϕ̃ + β) Ni dΓ (4.34)

Γ

∂n Γ

X Z Z

= ϕ̃ j σαN j Ni dΓ + σβ Ni dΓ . (4.35)

j Γ Γ

44

The evaluation of equation 4.35 for all i is equivalent to

Bϕ̃ + b, (4.36)

with a Matrix B and vector b. In principle, B is a mass matrix. Its dimension is not

that of domain Ω but of its surface Γ , though.

If the boundary conditions are included into the Θ-method scheme described by

equation 4.28, then the linear system of equations becomes:

K" n+1 + θ ∆t Kσ n+1 − Bn+1 ϕ̃ n+1 = (K" n − (1 − θ )∆t (Kσ n − Bn )) ϕ̃ n

+ ∆t θ bn+1 + (1 − θ ) bn . (4.37)

In the last equation the matrices B(·) and vectors b(·) are global, however, whereas

the previous equations describe an evaluation at the local level.

In Sec.3.2.1 the relevant differential equation for heat conduction problems has

been presented in the form:

∂ϑ

cp % − div (λ grad ϑ) = ẇgen . (4.38)

∂t

Thus, the weighted-residual Galerkin method leads to equation:

Z Z

∂ϑ

cp % − div (λ grad ϑ) w i dΩ = ẇgen w i dΩ. (4.39)

Ω

∂t Ω

Z Z Z

∂ϑ

cp % w i dΩ − div (λ grad ϑ) w i dΩ = ẇgen w i dΩ . (4.40)

∂t

|Ω {z } |Ω {z } | Ω {z }

1 2 3

The first integral 1 becomes:

P

Z Z

∂ϑ ∂ j ϑ̃ j N j

cp % w i dΩ ≈ cp % w i dΩ (4.41)

Ω

∂t Ω

∂t

X ∂ ϑ̃ j Z

= c p %N j w i dΩ. (4.42)

j

∂t Ω

45

The evaluation of this equation for all i can be written as a matrix-vector product:

∂~

ϑ

Mc p % . (4.43)

∂t

The heat capacity or mass matrix Mcp % is composed of the coefficients from evaluat-

ing the integral for the corresponding combination of i and j . If the Bubnov-Galerkin

method with w i = Ni is chosen, the mass matrix is symmetric, as index i and j can

be interchanged freely.

The second integral 2 is similar to the integrals in Sec. 4.2.1. As before, a

stiffness matrix is obtained, which will be called Kλ . Depending on the boundary

conditions, we may also get a matrix B and vector b. The evaluation of the third

integral 3 with heat sources ẇgen leads to a vector, which we will call wV .

Thus, the heat conduction equation can be represented algebraically as:

∂ ϑ̃

Mc p % + Kλ ϑ̃ − Bϑ̃ = wV + b. (4.44)

∂t

If the θ -method is chosen for time integration, the following linear system is

solved:

n+1

Mcp % n+1 + ∆tθ Kλ n+1 − Bn+1 ~ϑ = Mcp % n − ∆t (1 − θ ) (Kλ n − Bn ) ~ϑn

+ ∆tθ wV n+1 + bn+1

+ ∆t (1 − θ ) (wV n + bn ) . (4.45)

been restricted to Dirichlet, Neumann and Robin boundary condition. These condi-

tions can also be applied to problems of heat conduction. The Dirichlet boundary

condition is equivalent to imposing the temperature on some parts of the domain.

The other boundary conditions are applied through matrix Bλ and vector b. The

surface integral obtained together with Kλ is equivalent to a weighted heat flux

through the surface:

Z Z Z

n · (λ grad ϑ) w i dΓ = n · q̇ w i dΓ = q̇n w i dΓ . (4.46)

Γ Γ Γ

46

Hence, the Neumann boundary condition is equivalent to the imposition of a normal

heat flux density q̇n . It results in a vector b with component bi :

Z

bi = q̇n w i dΓ . (4.47)

Γ

The convection and radiation boundary conditions impose a normal heat flux

density, too. However, the flux depends on the temperature at the boundary and

is not known beforehand. Therefore, the convection boundary condition and the

linearized radiation boundary condition can be considered as special cases of the

Robin boundary condition.

In the case of convection, one has q̇n = h (ϑ − ϑ∞ ):

Z Z

q̇n w i dΓ = h ϑ̃ − ϑ∞ w i dΓ (4.48)

Γ Γ

X Z Z

= ϑ̃ j αN j w i dΓ + (−αϑ∞ )w i dΓ (4.49)

j Γ Γ

If this integral is evaluated for all w i and N j , one obtains a linear expression that

enters Bλ ϑ̃ + b.

The radiative heat flux boundary condition is imposed similarly:

Z Z

q̇n w i dΓ = εσSB ϑ̃4 − ϑsur

4

w i dΓ (4.50)

Γ Γ

Z

3 4

4

≈ εσSB (4 ϑ̃(−1) ϑ̃ − 3 ϑ̃(−1) − ϑsur w i dΓ (4.51)

Γ

Z Z

3

4

4

≈ 4εσSB ϑ̃(−1) ϑ̃w i dΓ − εσSB 3 ϑ̃(−1) + ϑsur w i dΓ . (4.52)

Γ Γ

which determines the irradiated heat. After linearization in the second line, the

radiation boundary condition has the form of a Robin boundary condition that

contributes to matrix Bλ and vector b The temperature ϑ̃(−1) designates the best

known value for ϑ̃, i.e., the solution of the previous iteration, when available.

47

5 Characterization of Nonlinear

Materials

5.1 Nonlinear Capacitances

In most materials electric field strength and field density are proportional, as

dielectric saturation does not occur before the breakdown of the dielectric [70].

Whenever these two fields are non-proportional, dielectric permittivity and electrical

capacitance become nonlinear. Varactor diodes are the classic example of a nonlinear

capacitor. They are exploited as such to tune resonant circuits. The width of their

depletion layer, which determines the capacity, depends on the applied voltage [66].

Another example are ferroelectric materials, which are not only nonlinear but show

in addition hysteresis effects.

Often one looks at the relationship between electric charge and voltage, their

q-u-characteristic q = f (u) [93]. For linear materials f (u) is simply described by

capacitance C as proportionality factor, i.e.:

q = C u. (5.1)

In analogy, one can use a voltage-dependent static capacitance C(u), such that

The electric current i is equal to the time derivative of the charge q̇:

dq du

i = q̇ = = Cd (u) u̇. (5.3)

du d t

Cd can be called differential capacitance. It describes the increase of charge for an

infinitesimally small voltage increase. When materials are linear, it is equal to the

static capacitance C . For nonlinear materials, C and Cd are distinct:

d dC ∂C

i= (C(u) u) = u + C u̇ = u + C u̇. (5.4)

dt dt ∂u

| {z }

Cd

48

guard ring

upper electrode

V

sample

lower electrode

Figure 5.1: Schematic diagram of the measurement setup. The system is excited by

an AC voltage signal. A guard ring surrounds the upper electrode to

eliminate the influence through stray fields on the current measured by

the ammeter. In consequence, the electric field is almost homogeneous

and the material parameters can be calculated, when the surface area of

the electrode and the thickness of the sample are known.

These definitions are not shared by all authors. For example, the above-cited

book by Lerch [66], prefers to define the nonlinear capacitance C in such a way that

the relationship i = C u̇ holds. For this reason the differential capacitance Cd of Eq.

5.3 is sometimes called nonlinear capacitance. One can also find the designation

small-signal capacitance.

In the context of this thesis, the definition of Eq. 5.2 is preferred. In particular, it

seems to be the more natural choice, because it maintains the connection between

MAXWELL’s equations and circuit elements for the nonlinear case. This will be shown

in the following section (Sec.5.2).

The materials are characterized using voltage and current data from a measurement

setup that is shown schematically in Fig. 5.1. The material sample, for example, a

disk of silicone rubber filled with microvaristors (in Sec. 5.4.3), is placed between

two cylindrical electrodes, which are connected to an a.c. voltage source. The upper

electrode is divided into an inner electrode and a surrounding guard ring. While

the current through the guard ring is affected by fringe effects, the field between

the inner part of the upper electrode and the lower electrode, which determines

the current measured with the ammeter, is almost homogeneous and perpendicular

49

upper electrode

guard ring

sample

lower electrode

Figure 5.2: Photograph of the most relevant part of the described measurement

setup. It shows the electrodes and a silicone rubber sample used for mi-

crovaristor measurements by M. Draude at the High-Voltage Laboratory.

Printed with permission.

50

to the electrode surfaces. The guard ring configuration, which has been used for

the characterization of microvaristor samples at the High Voltage Laboratory of TU

Darmstadt, is shown in the photograph of Fig. 5.2.

Inside the sample, the current between the two electrodes flows in almost vertical

direction. The total current density J is the sum of resistive and capacitive current

density Jr and Jc :

J = Jr + Jc (5.5)

∂

= σ(E)E + "r (E)"0 E . (5.6)

∂t

In this equation, both conductivity σ and relative permittivity "r are allowed to

be nonlinear and depend on present field strength. The scalar field strengths and

current densities correspond to the non-zero component of the field strength vector.

Since the field is almost homogeneous, the voltage difference u is equal to u = Ed

with thickness of the sample d . The measured current i corresponds to the relevant

surface area A (of the upper electrode in Fig. 5.1) times current density J , i.e.,

i = JA.

Based on these assumptions one can derive the voltage-dependent conductance

and capacitance:

i = (Jr + Jc )A (5.7)

u ∂ u

= σ(u/d) + "r (u/d) "0 A (5.8)

d ∂t d

A ∂ A

= σ(u/d) u + "r (u/d) "0 u (5.9)

d ∂t d

d

= G(u) u + (C(u) u) . (5.10)

dt

Thus, the capacitance C(u) is equal to the capacitance of a parallel-plate capacitor

with surface area A and distance d between the plates. The only particularity is that

the relative permittivity "r is allowed to depend on electric field strength:

A

C(u) = "r (u/d)"0 . (5.11)

d

Similarly, the conductance G(u) is that of a straight wire with uniform cross-section

area A and length d :

A

G(u) = σ(u/d) . (5.12)

d

(5.13)

51

Inversely, if the above equations hold for G(u) and C(u), then:

d

"r (E) = C(Ed) , (5.14)

"0 A

d

σ(E) = G(Ed) . (5.15)

A

These two equations show how the material properties are related to the proper-

ties of the bulk material that are measured.

from measurements of voltage ui and current ii . A measurement setup like the

above-described is excited by an a.c. voltage signal. Both voltage and current are

recorded. Given that the geometry corresponds to that of Sec. 5.2, the relationships

of Eqs. 5.11 to 5.15 hold.

10 du 500

dt

=0

ir

5 250

Voltage / kV

u=0 Current / µA

0 0

−5 ic −250

−10 −500

−20 −10 0 10 20

Time / ms

Figure 5.3: Illustration for the method of chosen points. At zeros and local extrema

of the voltage signal, the current can be considered either completely

resistive (i = ir ) or capacitive (i = ir ). The shown electric field strength

(blue) and current density (red) are derived from the measured voltage

and current of a silicone rubber specimen filled with microvaristors at

high peak voltage.

52

For a linear material, the unknown values for linear conductance G and capaci-

tance C can be derived from two data points P1 = (u1 , u̇1 , i1 ) and P2 = (u2 , u̇2 , i2 ) by

solving the linear system of equations

i1 = Gu1 + C u̇1 (5.16)

i2 = Gu2 + C u̇2 , (5.17)

with measured voltage ui and current ii .

This kind of approximation can be called method of chosen points [92]. The

choice of points, for which either ui = 0 or u̇i = 0 is fulfilled, makes the evaluation

particularly simple. When voltage reaches its peak value, the entire current is

resistive, implying G = i/u. At voltage zero-crossings, the current is exclusively

capacitive, such that C = i/u̇.

Although in practice more points might be used to reduce noise, it is in principle

possible to obtain the material parameters simply by evaluating

G = i1 /u1 with u̇1 = 0 (5.18)

C = i2 /u̇2 with u2 = 0. (5.19)

Fig. 5.3 illustrates this approach which can still be applied in the case of a

voltage-dependent electrical conductance, even though one searches a number of

conductance values G(ui ) that are valid for n different voltages ui instead of one

constant value G :

G(ui ) = ii /ui with u̇i = 0 (5.20)

C = in+1 /u̇n+1 with un+1 = 0. (5.21)

However, the application of this method to materials with equally nonlinear

permittivity leads to erroneous results like those discussed in Sec. 5.4.3. The

error between observed and predicted currents due to the nonlinearity of relative

permittivity amounts to:

∂C

ierr = C(u) − C(0) + u u̇. (5.22)

∂u

The interesting point about this error is that it does not only depend on the excitation

voltage, but it matters, if voltage is rising or falling. If one plots the difference

between measured and calculated current as function of voltage, on may observe

loops suggesting the presence of a hysteresis phenomenon, even if that is not the

case.

Instead of extending the described method to materials with nonlinear conduc-

tivity and permittivity, another, more general approach was developed that will be

presented in the following section (Sec. 5.4).

53

5.4 Least-Squares Based Method

As alternative to the method of Sec. 5.3, another approach has been tried and

resulted in more than satisfactory results for the characterization of materials

with nonlinear permittivity and conductivity (see Sec. 5.4.3). In this alternative

approach the material parameters are not obtained by comparison of data from

individual sampling points which fulfill specific conditions. The data are instead used

collectively to obtain nonlinear material characteristics that minimize statistically the

error between observed and predicted currents. Although the approach proves its

usefulness and is relatively trivial, the approach seems not to have been used earlier

in electrical engineering. It can be compared to the identification of spatially-varying

material properties that is used for many years in soil hydrodynamics, e.g., [60].

The method seeks to obtain a solution that directly satisfies:

d

i = G(u) u + (C(u) u) . (5.23)

dt

∂

J = σ(E)E + ("r (E)"0 E) . (5.24)

∂t

coefficients σ̃k and their corresponding basis functions wσk (E):

X

σ ≈ σ̃ = σ̃k wσk (E). (5.25)

k

Analogously:

X

"r ≈ "˜ = "˜k w"k (E). (5.26)

k

X

J̃ r = σ̃k wσk (E) E. (5.27)

k

54

This can be rewritten in the form of a matrix-vector product. If the field strength

was measured for n sampling moments and k coefficients are used to describe σ̃i ,

then:

r wσ (E )E

wσ2 (E1 )E1 · · · wσk (E1 )E1 σ̃

J̃1 1 1 1

1

..

J̃2r

w σ

(E )E w σ

(E )E .

σ̃2

. = 1 2 2 2 2 2

. .. .. .. .. . (5.28)

. . . . .

J̃nr wσ (En )En ··· · · · wσk (En )En σ̃k

| 1 {z

}

G

∂

J̃ c = ("r "0 E) (5.29)

∂t

∂ X

= "˜k w k (E) "0 E

"

(5.30)

∂t k

X

∂ w"k (x) ∂E

= "˜k w k (E) +

" ·E . (5.31)

k

∂ x x=E ∂t

The resulting linear system stretches over more than the width of a column. Thus,

it is not displayed in the same way as Eq. 5.28. Nevertheless, a comparable linear

system is obtained:

J̃c = C "˜r . (5.32)

σ̃

J̃ = A , (5.33)

"˜

where A is a matrix composed of the submatrices G and C , i.e., A = G|C . The

subvectors σ̃ and "˜ contain the unknown coefficients, which describe the material.

The matrix A is asymmetric and possesses generally much more rows than

columns. While the number of rows corresponds to the potentially very large

number of retained sampling points, the number of columns is equal to the number of

coefficients deemed necessary to characterize conductivity and relative permittivity.

Furthermore, the system may be sparse, depending on the choice of basis functions.

The goal is getting a set of coefficients for σ̃ and "˜, for which the difference

between observed and estimated current is minimized:

55

minimize
J̃(σ̃, "˜) − Jobs
(5.34)

σ̃,˜

"

σ̃

AT A = AT J, (5.35)

"˜

which provides an ordinary least squares (OLS) estimate of the coefficients [36, 98].

In practice, one often prefers not to solve this equation and uses a singular value

decomposition (SVD) or QR decomposition of matrix A, because of the increased

condition number of the normal equation.

5.4.2 Implementation

tion that the material parameters do not depend on the direction of the electric field,

but only on its absolute value. Furthermore, non-uniform B-splines, or basis splines,

of arbitrary order were chosen as basis functions. Following de Boor [14, 15], the

B-spline of order 1 for a sequence of knots x := (x i ) is defined by

¨

1, if x i ≤ x < x i+1 or x = max j (x j )

Bi1 (x) := (5.36)

0, otherwise.

x − xi x i+k−1 − x

Bik := Bi(k−1) + B(i+1)(k−1) . (5.37)

x i+k−1 − x i x i+k−1 − x i

56

Bi,1 Bi+1,1 Bi+2,1

1

0

xi x i+1 x i+2 x i+3

Bi−1,2 Bi,2 Bi+1,2 Bi+2,2

1

0

xi x i+1 x i+2 x i+3

1

Bi−1,3 Bi,3 Bi+1,3

0

1 xi x i+1 x i+2 x i+3

Figure 5.4: The diagram shows B-splines Bik of orders one to three, according to the

definition by de Boor, for the same n knots x i . The upper part shows

B-splines of order 1, which are piecewise constant. In the middle part,

one sees B-splines of order 2, which are also known as hat functions or

linear B-splines. The lower part shows B-Splines of order 2, or quadratic

B-splines.

Some splines of different degree are shown in Fig. 5.4. If the sequence of knots

is allowed to contain some knots multiple times, this special case has to be treated

differently. However, we will require that the nodes are in strictly ascending order

Multiple knots are equivalent with reduced continuity requirements at these knots,

but for the numerical field simulations continuity of the material characteristics at

least up to their first derivative is desired, even if the material possessed indeed

some discontinuity. In this case, one obtains a spline curve s(x) as weighted sum of

the B-splines with weights ai :

X

s(x) = ai Bik (x). (5.38)

i

57

This curve is piecewise-polynomial of degree k − 1 between two neighboring

knots. At the knots themselves, the curve is continuous up to its (k − 2)-th derivative.

In our context, s(x) corresponds to either σ̃(|E|) or "˜r (|E|). For the purpose of

estimating the material parameters, the set of knots spans the interval between

0 V/m and the maximum observed field strength Ê . This interval is subdivided into

n smaller intervals, which form a set of (n + 1) knots [33]:

Not requiring a uniform distribution, the knots can be chosen freely as long as

they are sorted in ascending order. By default the nodes are placed at equal distance

between minimum and maximum value. During the calculation of the B-splines,

some additional knots are adjoined temporarily, to remove the imposition of conti-

nuity at the upper and lower limits of the estimation interval, which unnecessarily

restricts the shape of the material characteristics.

Using B-splines of order k as basis functions, some set of knots ∆ and materi-

als that depend on the present absolute value of field strength, the resistive and

capacitive current densities become

X

J̃ r = σ̃i Bik

σ

(|E|) E (5.40)

i

and

X

∂ Bik

"

(x) ∂E

c

J̃ = "˜k Bik

"

(|E|) + |E| . (5.41)

i

∂ x x=|E| ∂t

The two equations are the equivalent of Equations 5.27 and 5.31 and the rest of

the procedure is the same as described above.

The time-derivative of field strength ∂∂ Et is of critical importance. Calculation

by means of finite differences from measured time and voltage leads to large

oscillations. The application of moving averages and low-pass filters reduces the

problem. However, the best results were obtained by calculating an approximating

polynomial spline curve and evaluating its analytical derivative at the sampling

points [42, 25]. The estimation of this spline curves is done in essentially the same

way as the subsequent estimation of the material parameters.

silicone rubber material. S. Blatt, who had realized the measurements and kindly

58

provided his data, observed that the characteristics, that he obtains with the method

of Sec. 5.3, deliver an unsatisfactory explanation of the material behavior and ques-

tioned the linearity of permittivity. Others have suggested a nonlinear permittivity

of varistor materials, too [20], but without providing strong enough arguments.

As mentioned above, the conventional method, according to Sec. 5.3 faces diffi-

culties for materials with nonlinear conductivity and permittivity. These difficulties

have been the motivation for the present approach. In the following, an OLS-based

method will be applied to estimate conductivity and permittivity. The resulting ma-

terial characteristics offer a much better agreement between observed and expected

current. Furthermore, they confirm the argument for a nonlinear permittivity of

varistor materials.

Fig. 5.5 illustrates the deficiency of the model with field-dependent conductivity

and constant permittivity. Normally, the measured current should closely resemble

10 500

5 250

Voltage / kV

Current / µA

0 0

ir 2

−5 −250

ir 1

u1 = u2

−10 −500

−20 −10 0 10 20

Time / ms

Figure 5.5: This figure shows that the assumption of a constant permittivity and a

field-dependent conductivity can not be upheld for the microvaristor-

filled specimen. The currents and the slope of the voltage signal at the

voltage zero crossings in Figure 5.3 lead to a relative permittivity of about

"r = 18.35. The corresponding capacitive current is subtracted from

the observed current to obtain the resistive current (dashed line). If

the model was valid, the same currents should be observed for identical

voltages. Notwithstanding u1 = u2 , the current i r1 is much larger than

ir2 .

59

60

Conductivity / nS/m

50

40

30

20

10

0

0 5 10 15

Field Strength / kV/cm

23

22

Rel. Permittivity

21

20 ǫr = const.

19

18

17 ǫr = ǫr (|E|)

16

0 5 10 15

Field Strength / kV/cm

Figure 5.6: Estimated electrical conductivity σ and relative permittivity "r of a mi-

crovaristor sample. The characteristics have been obtained using the OLS

approach for two different cases: a constant relative permittivity (violet),

and a field-dependent relative permittivity (blue). The estimated conduc-

tivity seems to be independent from the assumption on permittivity, so

that the violet and blue characteristics are almost identical. Conductivity

increases sixfold. By contrast, permittivity is shown to grow by almost

50 % as field strength increases and differs clearly from the curve for a

relative permittivity.

60

Peak Voltage: 1.5 kV Peak Voltage: 6 kV

17.3 69.1

10 20 30 40 50 60 10 20 30 40 50 60

−17.4 −69.4

Current Density / mA/m2

34.3 102

10 20 30 40 50 60 10 20 30 40 50 60

−39.8 −101

Peak Voltage: 4.5 kV Peak Voltage: 9 kV

52.5 172

10 20 30 40 50 60

10 20 30 40 50 60

−74 −172

Time / ms

Figure 5.7: The diagrams show measured (black) and predicted current densities

according to the characteristics of Fig. 5.6. The blue curves represent the

predicted current densities for a constant relative permittivity "r , whereas

the predicted current densities for a nonlinear relative permittivity are

drawn in red. The assumption of a nonlinear permittivity leads to a much

better agreement between observed and predicted behavior.

61

the sum of resistive and capacitive currents. After subtraction of the estimated

capacitive current, the residual current has completely different values for the

same voltage. This is contradictory to an interpretation of the remaining part as

resistive current. Therefore, the model with constant permittivity and field-strength

dependent conductivity poorly explains the observations.

The measured data, for which results are presented, comprise time, voltage and

current for a series of thirty-five measurements with different peak voltage. The

setup corresponds to that described in Sec. 5.2. The peak a.c. voltage of the

individual measurements varies in steps of 250 V between 500 V and 9000 V . Since

the examined specimen has a thickness of about 5 mm, the maximum field strength

is 18 kV/cm. For every voltage level, the data consist of 2000 sampling points that

span over two periods.1

The upper part of Figure 5.6 shows the characteristic curves for electrical conduc-

tivity, while the lower part of the figure depicts the corresponding curves for relative

permittivity. The two parts of the figure plot the characteristics under the alternative

assumptions of a constant relative permittivity (in violet) and of a field-strength

dependent permittivity (in blue). The blue and violet curve are almost identical in

the upper part of the figure. Apparently the estimated conductivity is independent

from the assumptions made about relative permittivity. In both cases, conductivity

increases approximately sixfold with rising field strength.

The lower part of the figure shows relative permittivity. Requiring a constant

relative permittivity, its value is estimated to be close to 20. If permittivity is

allowed to vary with field strength, the blue curve is obtained. It shows that relative

permittivity increases by almost 50 %.

In Fig. 5.7 it is shown how good the different models fare in explaining the

observed data. Given the evolution of field strength, derived from the measured

voltage difference between the two electrodes, and the estimated characteristics for

conductivity and relative permittivity, one can calculate the respective predicted cur-

rent densities. Both models reproduce up to a certain point the observed current or

current density. However, the difference between measured and predicted behavior

is much larger for the model with constant relative permittivity. The assumption

of a nonlinear permittivity allows a significantly better prediction of the observed

currents. This is noticeable in the reduced overestimation of peak current density

for the measurements with low peak voltage and the better approximation of the

signal shape at very high field strengths.

1

Actually, the input data are already post-processed data from measurement data over a larger

time span.

62

Peak Voltage: 1.5 kV Peak Voltage: 6 kV

24.0

29.0

18.0

17.1

10 20 30 40 50 60 10 20 30 40 50 60

Peak Voltage: 3 kV Peak Voltage: 7.5 kV

Error / % of Jmax

24.9

34.3

14.5

23.7

10 20 30 40 50 60 10 20 30 40 50 60

Peak Voltage: 4.5 kV Peak Voltage: 9 kV

38.7

31.6

26.4

15.6

10 20 30 40 50 60 10 20 30 40 50 60

Time / ms

Figure 5.8: Difference between calculated and measured current density (see

Fig. 5.7) with respect to maximum

measured current density, i.e.,

J̃ (t) − J̃ (t) / max J obs (t) . The error of the model with constant

i i t i

permittivity is drawn in blue, while the significantly lower error of the

model with nonlinear permittivity is shown in red.

In Fig. 5.8, the difference between predicted and actually observed current density

is presented as a percentage value with respect to the respective maximum observed

current density, which has been chosen as a reference to avoid spikes close to the

63

zero crossings. This difference or error is calculated according to the following

equation from the data already shown in the last figure:

J̃ (t) − J obs (t)

i i

errori (t) = . (5.42)

max J obs (t)

i

t

The plots show that the prediction error for the current density has been reduced

significantly for the measurements with low and high peak voltages. For the mea-

surement with 1.5 kV peak voltage, the predicted current density differs by up to

29 % of its maximum observed value. This value decreases to about 17.1 %, which

constitutes an improvement of more than 40 %. The prediction quality increases

by an even larger value, almost 60 %, for the measurement with the highest peak

voltage. This is particularly important, as understanding the behavior at high field

strengths is critical to manage electrical losses and the dissipation of heat. In the

middle range of peak voltages the improvement of the prediction quality is smaller,

but still substantial.

1

10−1

10−2

Order k=1

10−3 k=2

1

Rel. Error

10−4

2

−5

10 3

10−6 4 k=3

k=4

5 k=5

10−7

6 k=6

10−8

1 2 3 5 7 10 15 20 30

Degrees of Freedom

different order for increasing number of degrees of freedom. In all cases

the material model converges to the same relative error of about 12.07 %

according to Eq. 5.43. The value obtained for 30 B-splines of order six

constitutes the reference solution.

64

60

Conductivity / nS/m

50

40 f

30 e

20 d

c

10 a b

0

0 5 10 15

Field Strength / kV/cm

26

26

Rel. Permittivity

22

20 f

e

18 d

16 a b c

0 5 10 15

Field Strength / kV/cm

Figure 5.10: Characteristics for conductivity σ and relative permittivity "r of a mi-

crovaristor sample. The characteristics were obtained from measure-

ments with different AC peak voltage (a: 1500 V, b: 3000 V, c : 4500 V,

d : 6000 V, e: 7500 V, f : 9000 V). For comparison, the nonlinear char-

acteristics of Fig. 5.6 have been added as dashed line. Except at very

high field strength, the conductivity seems to depend more on peak

voltage than on the present field value, which suggests the presence of

a relaxation process.

65

Finally, Fig. 5.9 shows the convergence of the estimates with increasing number

of degrees of freedom for different B-spline orders. While Fig. 5.8 has already given

us some means to judge the error of the individual measurements, one can define

another global measure of the total error. For a constant permittivity the relative

error η over all observations, given by

J − J̃

obs

η= , (5.43)

kJobs k

converges towards 23.53 % with increasing number of degrees of freedom and spline

order. Accepting a nonlinear permittivity the relative error decreases quickly towards

a final value η0 ≈ 12.07%. The convergence diagram displays the convergence to

this value:

|η(d, k) − η0 |

, (5.44)

η0

number of degrees of freedom for both conductivity and permittivity, the total

number of estimated parameters is twice that value. Not knowing the precise value

of η0 , the supposedly most accurate estimate η(30, 6) is used as reference value.

This leads to the convergence behavior shown in the figure and is in good agreement

with the expected behavior, which is indicated by the dashed lines.

The results prove that the least-squares based method is applicable to estimate the

electrical properties of nonlinear materials and provides good results. The obtained

characteristics indicate strongly that relative permittivity is indeed nonlinear and

that this effect plays an important role in explaining the observable behavior during

measurements.

These significant results have been the point of departure for further studies after

realizing that the material is not only nonlinear, but that there is also some kind of

memory effect involved. One should expect that the residual difference between

predicted and observed values is random and that the exclusion of individual

measurements does not affect the estimated characteristics in the covered field-

strength range. However, this is not the case. It is shown in Fig. 5.10 that the

obtained characteristics differ significantly with respect to the retained measurement

data. For example, if one uses the data from the measurement with peak voltage

9000 V, one estimates characteristic curves for conductivity and relative permittivity

that lie substantially higher than those obtained from the other measurements.

If only the measurement data for the corresponding signal are used to estimate

the material characteristics, the resulting characteristic curves for conductivity and

permittivity explain the observed behavior consistently better than the previously

66

Peak Voltage: 1.5 kV Peak Voltage: 6 kV

13.7 67.3

10 20 30 40 50 60 10 20 30 40 50 60

−14.4 −67.9

Current Density / mA/m2

27.6 102

10 20 30 40 50 60 10 20 30 40 50 60

−32.1 −102

Peak Voltage: 4.5 kV Peak Voltage: 9 kV

44.3 172

10 20 30 40 50 60

10 20 30 40 50 60

−60.3 −176

Time / ms

Figure 5.11: Measured and predicted current density for the characteristics of Fig.

5.10. The application of different characteristics depending on the

respective peak Voltage leads to an important improvement of the

prediction quality. It has become difficult to spot the measured current

density (black) behind the current density calculated for the model

with peak-voltage dependent nonlinear permittivity and conductivity

(red). The blue curves show the predicted current density under the

assumption of a constant permittivity.

67

Peak Voltage: 1.5 kV Peak Voltage: 6 kV

4.0 15.2

3.0

3.4

10 20 30 40 50 60 10 20 30 40 50 60

Peak Voltage: 3 kV Peak Voltage: 7.5 kV

Error / % of Jmax

20.3

8.8

3.6

2.9

10 20 30 40 50 60 10 20 30 40 50 60

Peak Voltage: 4.5 kV Peak Voltage: 9 kV

12.9

29.8

8.5

5.7

10 20 30 40 50 60 10 20 30 40 50 60

Time / ms

Figure 5.12: The individual figures provide an illustration of the error between ob-

served and predicted currents in Fig. 5.11. The model with nonlinear

conductivity and permittivity (red) proves to provide a much better

prediction than the model with nonlinear conductivity and constant

permittivity. The values should be compared with the previous Fig. 5.8,

which gives the same representation, though for the case of a common

characteristic irrespective of peak voltage.

68

1334

1000 500

300

Conductivity / nS/m

200

100

100

50

10

30

20

10 15 10

10

5

0 0

Peak Field Strength / kV/cm Field Strength / kV/cm

26

26 24

24

Rel. Permittivity

22 22

20

18 20

16

18

10 15

10 16

5

0 0

Peak Field Strength / kV/cm Field Strength / kV/cm

silicone rubber as function of present and peak field strength, as derived

from a series of 50 Hz measurements. The conductivity is shown to

depend not only on present field strength, but also exponentially on

the observed peak field strength of the respective measurement.

69

200 Polynomial Order

æ

100 5

50 2

Rel. Error / %

20

æ

æ

1

10

5 æ

2

1

5 10 15 20 25 30 35

Number of Retained Measurements

The horizontal axis shows how many of the 35 measurements are used

for the characterization of the material, while the vertical axis gives an

error (see Eq. 5.43) for the whole group of measurements. The color of

the points indicates the polynomial degree in peak voltage direction for

the given data point.

presented characteristics. This is shown in Fig. 5.11. While the difference between

red and black curves has already been small in Fig. 5.7, the two curves are now

almost indiscernible. It has become difficult to spot the black curve behind the red

one.

As before, the error is shown more explicitly in a second figure. Fig. 5.12 shows

that the relative error with respect to the peak current density of the measurement

is significantly reduced by taking into account present field strength and peak

voltage. For example, one can observe that the error is reduced from up to 15.6 %

to values below 5.7 %, if one looks at the critically important behavior at highest

field strengths, i.e., the measurement data for 9000 V peak voltage.

It has already been observed that the adoption of a nonlinear relative permittiv-

ity leads to material characteristics which are in much better agreement with the

experimental observations than those obtained under the assumption of constant

permittivity. Nevertheless, the explanation of the material behavior is still incom-

plete, as the influence of peak voltage, respectively peak field strength, shows. The

70

observed material behavior can be explained even better, if peak voltage is taken

into account, too.

The presented method can be extended to include a second dimension. This had

been done to obtain field-strength and temperature dependent characteristics (see

Sec. 6.1), but it is equally possible to use peak field strength Ê as second variable.

This leads to:

XX

σ = σ(|E|, Ê) ≈ σ̃w i (|E|)v j ( Ê) (5.45)

i j

XX

"r = "r (|E|, Ê) ≈ "˜r w i (|E|)v j ( Ê), (5.46)

i j

with basis function w i and v j which are not necessarily the same for conductivity σ

and relative permittivity "r .

Some characteristics obtained in this way are shown as surface plots in Fig. 5.13

(for 15 × 12 intervals with polynomials of degree 5 × 5). This kind of characteris-

tics are expected to provide accurate simulation results, when the electric in the

examined microvaristor material oscillates indeed with a frequency of 50 Hz, the

frequency for which the characterization data were obtained.

That it is possible to interpolate the obtained material characteristics to calculate

currents and fields at further voltage or field strength levels is demonstrated by

Fig. 5.14. It shows the minimum error according to Eq. 5.43 achieved for a specified

maximum of 320 degrees of freedom, 16 in field strength direction times up to ten

in peak field strength direction for conductivity and permittivity. Using every fifth

measurement is sufficient to characterize the material with an error that is below

5 % according to the previously used measure of error. With increasing number of

measurements taken into account to characterize the material the error converges to

less than 2.9 %. Furthermore, the figure shows that the choice of splines of second

or third order, equivalent to segment-wise linear or quadratic polynomials with

respect to peak field strength, is preferable to higher-order basis functions for the

same number of degrees of freedom.

In many cases, when the field is indeed oscillating with 50 Hz, this model may be

useful. If one is interested in an analysis at another frequency, one needs another

pair of characteristics, though. And for arbitrary signals, such a model is simply

inadequate. Therefore, the next section discusses relaxation and the attempts to

find a generally valid material model.

5.5 Relaxation

In the preceding sections the least-squares based approach had been successfully

applied to characterize a microvaristor material. It had even been possible to

71

generate distinct characteristics for each of the 35 a.c. voltage measurements.

Unexpectedly, the curves are not similar. Apart from differing with respect to the

maximum covered field strength, they show a clear upward shift of conductivity

and, to a lesser extent, of permittivity with respect to the peak voltage of the

measurement (see again Fig. 5.11).

As an increase in temperature or any other measurement-related influence was

excluded by colleague S. Blatt, one had to conclude that some kind of relaxation

phenomenon is involved. The material properties depend not only on momentary

field, but in some unknown way on the past evolution of the field.

The presentation of the constitutive relationships for Maxwell’s equations in

Sec. 3.1.2 includes an introduction into relaxation and its connection with dispersion.

While relaxation can normally be deduced from an analysis in frequency domain,

this is not possible in the present case. As the material is nonlinear with respect to

field strength, the transformation between frequency and time domain is no longer

allowed.

Nevertheless, one may assume that the microvaristor relaxation behavior is not

completely different and try to use a similar model to analyze its behavior. The

relaxation model of Eq. 3.16:

Z∞

D(t) = "0 "∞ E(t) + Φ̇(τ)E(t − τ) dτ . (5.47)

0

the observed relaxation is due to permittivity or conductivity, the same approach is

used for permittivity and conductivity:

Z∞

D(t) = "0 "∞ (|E(t)|) E(t) + "s (|E(t − τ)|) e−τ/τ" E(t − τ) dτ (5.48)

0

Z ∞

J(t) = σ∞ (|E(t)|) E(t) + σs (|E(t − τ)|) e−τ/τσ E(t − τ) dτ. (5.49)

0

This approach corresponds to the Debye relaxation model, if the material proper-

ties are not affected by field strength. For this approach the estimated total current

density becomes:

Z∞

J̃(t) =σ∞ (|E(t)|) E(t) + σs (|E(t − τ)|) e−τ/τσ E(t − τ) dτ

0

Z ∞

∂

+ "0 "∞ (|E(t)|) E(t) + "s (|E(t − τ)|) e−τ/τ" E(t − τ) dτ . (5.50)

∂t 0

72

Model η

"r = const. 23.53 %

σ, "r = f (|E|) 12.07 %

σ, "r = f (|E|, Ê) 2.89 %

relaxation 9.74 %

Table 5.1: Comparison of the attained error level with several material models. The

error is significantly reduced by accepting a nonlinear permittivity and is

further reduced, when peak field strength is taken into account as second

parameter. For the described relaxation model, an error η of about 9.74 %

has been achieved, which is better than the 12.07 % of the model relying

only on present field strength, but still far from the 2.89 % of the model

with peak field strength as parameter.

The parameters, on which the estimated current density J̃ relies, contained in the

characteristics "∞ (|E(t)|), σ∞ (|E(t)|), "s (|E(t)|) and σs (|E(t)|) as well as time

constants τ" and τσ .

In principle, the material properties are estimated in the same way as in Sec. 5.4,

i.e., by minimization of the difference between measured and estimated current

densities:

minimize
Jobs − J̃("∞ , σ∞ , "s , ςs , τ" , τσ )
2 (5.51)

"∞ ,σ∞ ,"s ,σs ,τ" ,τσ

Although the problem looks similar, apart from the larger number of parameters,

it is much harder to solve. Because of the time constants which enter the model as a

multiplicative factor, the problem is no longer linear.

Instead of solving some linear system of equations to find the optimum solu-

tion, an iterative method has to be used to solve the multidimensional nonlinear

optimization problem. The Nelder-Mead simplex method, which was used, does

not guarantee the discovery of a global minimum, but may converge to local

minima [79].

In practice, the search for a optimum has been realized by gradually increasing

the number of unknowns and manual intervention once the algorithm seemed to

have reached a local minimum. This procedure does not guarantee finding the global

minimum and there is no reason to assume that it has been found for the examined

microvaristor material. For the calculation of the delayed material properties the

periodicity of the signal had to taken into account to evaluate the integrals, as

the expected relaxation times were not small with respect to the duration of the

measurement data.

73

14

12

σ∞ / nS/m

10

8

6

4

2

0

0 5 10 15

Field Strength / kV/cm

30

25

20

ǫ∞

15

10

5

0

0 5 10 15

Field Strength / kV/cm

120 000

100 000

σs / nS/(m s)

80 000

60 000

40 000

20 000

0

0 5 10 15

Field Strength / kV/cm

Figure 5.15: Estimated conductivity σ∞ , relative permittivity "∞ and the charac-

teristic for σs , which describes the time-delayed resistive current. The

characteristic for "s is not shown because of its minor contribution. The

time constant associated with σs is τσ ≈ 1.02 ms.

74

The Table 5.1 shows that the best obtained relaxation model achieves a reduction

of the error η (from Eq. 5.43) from 12.07 % to 9.74 %. This improvement is far from

negligible. However, the relaxation model explains the observed material behavior

much less than the model with peak field strength as variable. Unfortunately, the

cost in computation time and the tendency of the solution process to end up in local

minima do not let us firmly conclude, whether another model is needed or if the

model is valid and only the numerical solution unsatisfactory. Much better material

characteristics than those shown in Fig. 5.15 are likely to exist. For that purpose,

a better search method would be of great usefulness. The relaxed permittivity is

omitted in Fig. 5.15 because its contribution to the total current is negligible. With

time constant τσ ≈ 1.02 ms, the current corresponding to σs is mostly of greater

importance than the current related to σ∞ .

75

6 Surge Arrester Simulations

6.1 Computation of Stationary Operational Condition of Surge Arresters

surge arrester under operating conditions, when the arrester is connected to an

a.c. voltage source. The voltage may be either the normal operational voltage or

an elevated voltage used for testing the thermal stability of a surge arrester. The

Joulean heat

Heat dissipation

instability

threshold

ϑ↑

Power

ϑ↓

stable

operating point

unstable

stable

ϑ↑

Temperature ϑ

Figure 6.1: Illustration of thermal runaway. The diagram shows the Joulean heat

generation and heat dissipation of a valve element. The first intersection

point constitutes a thermally stable operating point. Up to the second

intersection point the element is thermally stable as dissipation of heat

to the environment exceeds heat input. Beyond the second intersection

point between the two curves the element is thermally unstable as the

injected heat is higher than the heat that is simultaneously dissipated to

the environment. Adapted from Lat [64].

76

purpose of this simulation is to determine the stationary temperature distribution of

the arrester and to verify that no thermal runaway takes place.

Varistors possess a so-called NTC behavior, which means that their electrical

resistance decreases for rising temperature. Since the arrester voltage is given, the

Joulean heat losses increase. Normally, this constitutes no problem. The arrester

will return to its equilibrium temperature distribution, because heat dissipation

increases faster than the Joulean losses. This holds true for some temperature range

around the thermal equilibrium. At higher temperatures the Joulean heat losses

begin to exceed heat dissipation, though, as illustrated in Fig. 6.1. If the feed-in

voltage is not reduced in time, the arrester continues to heat up until it fails.

Guaranteeing thermal stability is a problem that has gained importance in the last

few years along with the development of surge arresters for higher voltage levels.

The greater height of ultra-high voltage (UHV) surge arresters, which reach voltage

levels of 1100 kV to 1200 kV today, causes severe design problems because of the

greater stray capacitances [47].

A thermal runaway does not occur spontaneously, but happens after the injection

of a significant amount of thermal energy, e.g., by a series of lightning pulses.

According to Cooray [54] a typical ground flash consists of four to five lightning

strokes. First, a so-called stepped leader connects to the ground and initiates a

return stroke that travels from ground to the cloud base. After the return stroke

the current may continue to flow for some milliseconds. After the disappearance

of the current, the discharge provoked by the return stroke can evolve into a dart

leader traveling to ground through the previously established channel, resulting in

a further return stroke. More dart leaders and return strokes may follow.

Thus, one is concerned about the likelihood of events, which lead to thermal

runaway, and about the energy an arrester can absorb and still return to thermal

equilibrium.

This problem is not satisfactorily solved, yet. Zheng et al. [105] claim to have

succeeded in computing the thermal stability of a full-scale arrester, but rely on

major model simplifications, e.g., reducing the electrical part of the problem to

one dimension. The following section, Sec. 6.2, presents another approach with

its own limitations. A further, more conventional method is currently examined by

colleagues at TU Darmstadt.

If the problem of simulating the cooling behavior of a surge arrester over time

spans of several hours continues to be a problem, the determination of the thermal

equilibrium state seems to be more feasible. The idea is that of combining a transient

simulation of the electric part of the problem with a static solution of the thermal

problem:

− div (λ grad ϑ) = ẇgen . (6.1)

77

Start

ϕ0 (t = 0), ϑ0

ϕi (t = 0), ϑi

ẇgen,i

ϕi+1 (T )

ϑi+1

no

Converged?

yes

End

Figure 6.2: Flow chart illustrating the procedure to calculate a thermally stationary

state for a nonlinear, time-periodic electric problem, for which temper-

ature can be considered as constant. The time-averaged Joulean heat

losses are calculated by means of a transient electro-quasistatic simula-

tion. They are then used to get a temperature field by solving the static

heat conduction equation. To obtain a thermal and electrical steady

state after several iterations, the simulation is restarted with the thermal

and last electric solution.

78

10−5

10−5

Conductivity / S/m

10−5.5

10−6

10−6

10−6.5

10−7

400 2

1.5

1

0.5 10−7

300

Temperature / K Field Strength / kV/cm

tivity of a ZnO varistor.

The electric part is simulated for (at least) one period. During this simulation

the Joulean heat losses per period are calculated. Using these heat losses, the static

temperature distribution in the arrester can be determined. Since the losses in the

arrester depend on temperature, the electric simulation is repeated for the new

temperature. In parallel, the final distribution of the electric potential is used as

initial solution for the next iteration. In the course of several cycles the thermal

solution is expected to converge to its equilibrium distribution, while the electric

initial solution converges to its steady state (see Fig. 6.2).

Using this procedure the equilibrium temperature of a four-unit surge arrester

has been calculated, which is used regularly for experiments at the High-Voltage Lab

of TU Darmstadt. The simulation reproduces an experiment, for which the arrester

was placed without pedestal close to ground. This increases the nonhomogeneity of

the voltage stress and heat losses for a given voltage level. While this is undesirable

in practice, it can be useful for some experiments.

For conductivity and permittivity, field-strength and temperature dependent

characteristics are used, which have been obtained in the same way as the peak

field-strength dependent characteristics in Sec. 5.4.3. However, peak field strength

is replaced by temperature as second parameter in Eqs. 5.45 and 5.46.

79

1800

1600

Rel. Permittivity

1500 1400

1200

1000

1000

500

800

400 2 600

1.5

1

0.5 400

300

Temperature K Field Strength / kV/cm

of a ZnO varistor.

The data for the varistor characterization come from a large sequence of mea-

surements. They were provided by the High-Voltage Lab of TU Darmstadt. The

temperature spans over a range between about 20 and 200 ◦C. The characterization

process leads to estimated conductivity values in the range between 1 · 10−7 S/m

and 1 · 10−5 S/m, while the relative permittivity lies between 400 and close to 2000

(see Figs. 6.3 and 6.4).

To obtain better simulation results two boundary conditions are added. First, a

combined convection/radiation boundary for the surface of the arrester (housing,

flanges). The second boundary condition is of greater interest. It is a consequence

of a suggestion by M. Tuczek that radiative heat exchange between the resistor

column and the interior surface of the housing might be of importance. A complete

description of the interaction between all the elements on the surface of resistor

column and housing would lead to dense matrices, apart from posing questions

about the respective line of sight. Therefore, a simplified boundary condition for

the reciprocal heat exchange has been implemented to take into account the heat

exchange in radial direction.

The resistor column is subdivided into numerous segments, essentially the sur-

faces of the individual varistors. Each of these segments is supposed to exchange

80

356.4 K 488 kV

293 K 0 kV

Figure 6.5: Stationary heat distribution and momentary electric potential of a surge

arrester model. The model corresponds to a four-unit porcelain-housed

surge arrester used for experiments at the High-Voltage Lab of TU Darm-

stadt. Voltage stress and heat generation are increased , because the

arrester is not placed on a pedestal, as is normally the case. Boundary

conditions have been chosen in accordance with the guidelines for the

determination of the voltage distribution along metal-oxide surge ar-

resters of standard IEC 60099-4 [51]. The computational domain extends

in radial and upward direction beyond the shown area.

81

2.0

Field Strength / kV/cm

1.5

1.0

0.5

0.0

0 1 2 3 4

Distance from Arrester Base / m

Figure 6.6: Electric field strength along the resistor column for the arrester model

shown in Fig. 6.5. The voltage stress is greatest for the two topmost

arrester units. The curve seems to indicate that field grading succeeds in

homogenizing the voltage stress in the upper part of the arrester.

82

350

340

Temperature / K

330

320

310

300

290

0 1 2 3 4

Distance from Arrester Base / m

Figure 6.7: Temperature along the resistor column. At the bottom, the temperature

remains close to the assumed ambient temperature of 293 K, whereas

the equilibrium temperature exceeds 350 K in the upper arrester unit.

heat with a corresponding segment of the interior of the housing. Normally, radiative

heat exchange is modeled by:

q̇ = εσSB ϑ14 − ϑ24 , (6.2)

another surface at temperature T2 . The second temperature is generally given and

often corresponds to ambient temperature. In the present case, both temperatures

are unknown. The emission of heat can be calculated normally for the surface of

each finite element, whereas the irradiated heat is approximated from most recent

temperature solution of the surface with the corresponding boundary condition.

For each of the reciprocal heat surface segments, the total emitted heat power

is calculated by summing up the emitted heat power of the finite element surfaces

which belong to this segment. An average irradiation power density is calculated

by dividing this heat power through the surface area of the corresponding segment.

This is done before the loop through all elements to build the local system matrices.

83

During the loop the emission of heat is calculated element-wise as usual, while

the irradiation is imposed as additional Neumann boundary condition term. The

greatest deficiency of this model is that it ignores the impact of reflections and

assumes that the emitted heat from one surface is absorbed by the other one.

The simulations indicate that it is indeed possible to calculate the equilibrium

temperature in the presented way. Notwithstanding the nonlinear nature of the

problem, the same solution is obtained for largely varying initial conditions. If the

initial temperature is not chosen too high, so that the solution process converges

towards a thermal runaway solution, the same thermal equilibrium is attained. In

Fig. 6.5 the equilibrium temperature distribution is shown along with the electric

potential for the moment of peak excitation voltage. The temperature outside of the

arrester is meaningless, though, as convective heat is not passed from the arrester

to the air, but leaves the system at the housing and flange surface. The electric field

strength along the resistor column can be seen in Fig. 6.6 and the temperature along

the same line is shown in Fig. 6.7.

Regrettably, the results may not be compared directly with measurements. The

boundary condition, a Dirichlet boundary on the circumferential surface of the

simulation area, results in quite different results from those, which are obtained

in the real, three-dimensional environment. Furthermore, the results may vary

to a large extent with respect to variations of the thermal parameters, for which

uncertainty is great. Nevertheless, the values seem to be reasonable. They are also

in agreement with recent transient long-term stability simulations.

After the injection of thermal energy into a surge arrester by a lightning stroke or

a sequence of closely followed strokes, the arrester does not always return to its

previous stationary state. Sometimes a thermal runaway takes place. Unfortunately,

it is not immediately after the injection of thermal energy clear, whether the system

is still inside the limits of thermal stability. Even after waiting a long time, during

which the temperature seems to normalize, a thermal runaway may still materialize.

If one attempts to simulate this process, one faces the problem that one has to

follow the evolution of the arrester temperature over a long time span, possibly

many hours to safely exclude a thermal runaway. By contrast, the length of the time

steps, at least for the electrical solution, has to be about a millisecond to resolve

the exciting a.c. voltage signal and even less for the nonlinear currents. Such a

simulation is currently not feasible.

A possible escape from this difficulty is the application of an envelope equation

model. Such models are frequently used in optics and beam dynamics, e.g., [75, 35].

84

If one assumes that the electric potential inside the arrester is harmonic and sine

shaped, the potential can be substituted by a complex phasor Φ. Notwithstanding

that restriction, the amplitude and phase of the sine signal are allowed to vary

slowly in time and differ locally.

At every moment the momentary potential is given by

Φ(r, t) = Re Φ(r, t)e jωt . (6.3)

depend on location r and time t :

∂

div "r "0 grad Φ + div (σ + jω"r "0 ) grad Φ = 0. (6.5)

∂t

As long as Φ varies slowly, time steps many times as large as the period of the

harmonic signal are possible.

In the case of an arrester, one faces two questions with respect to this approach.

The first question refers to the assumption of an essentially harmonic potential. Due

to the nonlinearity of the varistors, it is not ensured that the real potential does

not deviate significantly from a sine shape. Even when such a signal is imposed at

its terminals, the potential is expected to deviate to some extent inside the resistor

column, because of stray capacitances and nonlinear varistor behavior.

Fortunately, a uniform distribution of the electric field along the resistor is a

major design goal for surge arresters. Therefore, one can hope that the deviation is

sufficiently small to justify the application of the envelope equation model.

The second question is related to the material parameters that one has to insert

into the equation. Every time step spans over many periods and during each of them

conductivity varies significantly along with field strength. The question is, which

value of conductivity has to be used to get a field solution that is close to the true

one. Certainly, the conductivity has to be in some way an average value of the real

conductivity in that time span.

Two different conductivities have been examined. The first conductivity relies on

a power equivalence: The Joulean losses per period with equivalent conductivity σP

are supposed to be equal to the losses with nonlinear conductivity σ(E):

Z T Z T

2

σP E (t) dt = σ(E(t))E 2 (t) dt, (6.6)

0 0

85

with period length T and sinusoidal field strength E(t) = Ê sin (ωt).

This leads to:

RT

0

σ(E(t))E 2 (t) dt

σP = RT . (6.7)

0

E 2 (t) dt

averaged current with a conductivity σI is required to be equal to the current with

real conductivity σ(E):

Z T /4 Z T /4

σI E(t) dt = σ(E(t))E(t) dt, (6.8)

0 0

the integration over a whole period leads to a value of zero.

Hence, the equivalent conductivity σI becomes:

R T /4

0

σ(E(t))E(t) dt

σI = R T /4 . (6.9)

0

E(t) dt

The potential of the envelope equation model is examined in two ways. For both,

the simplified model of multi-unit surge arrester from IEC standard 60099-4 is used.

The first test consists of a parameter study, where the local and total heat generation

of the envelope equation model is compared to the results of a normal time-domain

simulation. This comparison is realized for various material characteristics with

different degree of nonlinearity and voltage levels. The second study compares the

voltage stress along the resistor column with reference values cited in the mentioned

standard.

For the parameter study a varistor conductivity, which follows the power law kE α

plus an additional term σ0 for conductivity at low field strength is assumed:1

σ(E) = σ0 + kE α . (6.10)

1

Note that σ ∝ E α implies I ∝ U (α+1) , and not I ∝ U α .

86

α=0

3

Conductivity / µS/m

2.5

α=1

2

1.5 α = 2.5

α=5

1

α = 10

0.5

α = 30.89

0

0 0.5 1.0 1.5 2.0

Field Strength / kV/cm

the type σ(E) = σ0 +kE α . with a nonlinearity coefficient of α = 30.89 as

well as several fictitious characteristics with lower degree of nonlinearity,

but spanning the same range of conductivity values. These characteristics

were among those used to test the envelope equation model.

87

(σ̄I ,σ̄P )

25

Relative Error / %

20

(σ̄P ,σ̄P )

15

10

18 (σ̄I ,σ̄I )

17

16 5

15

14 0 (σ̄P ,σ̄I )

0 5 10 15 20 25 30

13

Nonlinearity Coefficient α

12 30

11

Power Density / kW/m3

10

(σ̄I ,σ̄P )

9

25

8

7

(σ̄P ,σ̄P )

20

6

5 (σ̄I ,σ̄I )

4

ref

3 15

(σ̄P ,σ̄I )

2

1

0 5 10 15 20 25 30

Nonlinearity Coefficient α

Figure 6.9: The upper diagram shows the relative error of the envelope equation

model for nonlinearity coefficients α between 0 and 30.89 for the simpli-

fied multi-unit surge arrester at measurement point 18 (left). The surge

arrester is excited with peak voltage Û = 570 kV, about 332 kV rms. The

names of the curves indicate, whether Eq. 6.11 or 6.12 has been used to

calculate the electric potential or Joulean heat generation. For example,

the solution obtained from using σI for electric potential and σP for heat

generation is named (σ̄I ,σ̄P ). The lower diagram shows the calculated

time-averaged heat losses for the different σI /σP -combinations as well

as the reference solution (ref).

88

This leads to equivalent conductivities σP and σI :

2kΓ 3+α

2

σP = σ0 + p Ê α , (6.11)

πΓ 2 + α2

p

πkΓ 1 + α2 α

σI = σ0 + Ê , (6.12)

2Γ 3+α

2

peak field strength Ê .

The fitting of the material characteristic of a real varistor, which had been realized

in a similar way to that described in Sec. 5.4, has led to the model:

30.89

E

σ(E) ≈ 2.07 · 10−7 + S/m. (6.13)

3.47 · 105 V/m

problems should manage to provide satisfactory results for nonlinearity coefficients

of α = 30 and beyond.

To examine the limits of the model, fictitious material characteristics were gener-

ated, which cover the same range of conductivity values, but exhibit a lower degree

of nonlinearity α (see Fig. 6.8). The multi-unit arrester was studied for various peak

voltages Û and conductivity characteristics. The time-averaged heat generation w̄

was calculated for a number of specified points in the resistor column, which are

shown on the left-hand side of Fig. 6.9. Normal electro-quasistatic simulations over

several periods to obtain a solution close to steady state provided reference values

for comparison with results for the envelope equation model.

In the simulations the peak voltage Û ranges from 235 kV to 564 kV and the

nonlinearity coefficient α from 0 to 30.89. For the highest simulated voltage levels

at 517 and 564 kV convergence and stability problems were observed for high values

of α leading to incomplete results. Therefore, the discussion will have to concentrate

on the results for a peak voltage of 470 kV, at which the problem is already markedly

nonlinear

Fig. 6.9 shows some of the simulation results for a peak voltage of 470 kV, which

corresponds to a rms voltage of about 332 kV. In the upper part of the figure the

relative difference of the calculated Joulean heat losses between the simulations with

the envelope equation model and the transient electro-quasistatic reference solution

can be seen. As expected, the results are practically identical for a nonlinearity

coefficient of α = 0. In that case the problem is completely linear and the envelope

equation model is exact, irrespective of the choice of equivalent conductivity.

89

18

17 423 kV 470 kV

16

15 282 kV

14 517 kV

18

13 376 kV

17 12 235 kV

Measurement Point

16

15 11 564 kV

14

13

10 329 kV

9

12

8

11

10 7

9

8 6

7

5

6 4

5

3

4

3 2

2

1 1

0 2 4 6 8 10 12 14

3

Power Density / kW/m

Figure 6.10: Volumetric heat generation due to Joulean losses at the 18 observation

points along the resistor column of the surge arrester shown on the

left. In an ideal arrester without the presence of stray capacitances the

electric field and heat losses are uniform. In practice, locally enhanced

fields are observed in the upper part of the arrester. The enhanced

electric fields, especially when they start exceeding the switching field

strength, lead to significantly higher losses, as can be observed for the

curves corresponding to peak voltages above 423 kV.

90

12 (σ̄P ,σ̄P ) 12 (σ̄P ,σ̄I )

470 kV

235 kV 470 kV

8 8

282 kV

376 kV 423 kV

Max. Rel. Error / %

4 423 kV 4

329 kV

0 0

0 5 10 15 20 25 30 0 5 10 15 20 25 30

12 470 kV

(σ̄I ,σ̄P ) 12 (σ̄I ,σ̄I )

8 8 470 kV

282 kV

4 4 235 kV 423 kV

376 kV

423 kV 329 kV 376 kV

0 0

0 5 10 15 20 25 30 0 5 10 15 20 25 30

Nonlinearity Coefficient α

Figure 6.11: Maximum relative error for each of the four equivalent conductivity

combinations as function of the nonlinearity coefficient α depending

on the labeled peak voltage. The maximum relative error is the highest

value among the 18 simultaneously observed relative errors. Generally

the maximum error is observed at a point in the upper arrester unit, but

not necessarily at the very top.

91

For α > 0 the error of the envelope equation model depends on the chosen

combination for the equivalent conductivity. The equivalent conductivity is used in

two contexts. First, it is used in the electro-quasistatic envelope equation model of

Eq. 6.5. Second, the conductivity is used to calculate the time-averaged Joulean

losses according to:

1

w̄ = σ̄i Ê 2 , (6.14)

2

where σ̄i stands for either σ̄I or σ̄P . In Fig. 6.9 the curves for relative error and

heat losses are designated according to the possible combinations. The first symbol

stands for the equivalent conductivity used for the field simulation and the second

one for the Joulean heat losses. For example, (σ̄I ,σ̄P ) stands for the result with

equivalent conductivity σ̄I used for the field simulation and σ̄P for heat generation.

According to the figure, the combination of using the power-equivalent conductiv-

ity for the field simulation and the current-equivalent conductivity for the Joulean

heat generation delivers the best results. The relative error is consistently below

10 % and is even lower than 3 % for high degrees of nonlinearity, The error of the

other approaches tends to rise with increasing degree of nonlinearity, although it is

still less than 10 % over the entire range of α-values for the combination (σ̄I ,σ̄I ).

It might have have been a reasonable guess to assume that the power-equivalent

conductivity is the preferable choice for the calculation of the power losses and that

current equivalence gives eventually a better field solution, but this is apparently

not the case.

Given that the field at the chosen observation point 18 close to the top of the

arrester is larger and more nonlinear than elsewhere, such relatively small errors

seem to indicate that the envelope equation model works reasonably well.

Below, the absolute value of volumetric Joulean heat generation is plotted for

the different equivalent conductivity combinations as well as the reference solution.

Unsurprisingly, the losses decrease with increasing degree of varistor nonlinearity.

Fig. 6.10 plots the volumetric Joulean heat generation of the reference solutions

at the 18 observation points, as shown on the left-hand side, for different peak

voltages between 235 and 564 kV. The previously discussed results correspond to

the topmost point of the 470 kV curve with losses exceeding 12 kW/m3 . Evidently,

this voltage is close to the protective level of the arrester. The resistance of the

resistor column decreases and current rises sharply.

The Fig. 6.11 offers a more complete visual representation of the errors of

the envelope equation model. It shows the maximum value of the relative error

observed along the resistor column for the four equivalent conductivity combinations

as function of nonlinearity coefficient α for peak voltages up to 470 kV. While

the maximum error is generally observed somewhere in the upper arrester unit,

92

18

2.0

12

Field Strength / kV/cm

1.5

6

1

0.5

0

80 85 90 95 100

Time / ms

0.2

Current Density / A/m2

0.15

18

0.1

0.05

12

6

0

80 85 90 95 100

Time / ms

Figure 6.12: Absolute value of field strength and resistive current density at various

points of the simulated arrester as function of time. The field-strength

data were obtained from the simulation for the reference solution at

Û = 470 kV with α = 30.89, while the current density values were

generated with Eq. 6.13. The curves correspond to measurement points

(6, 12 and 18) close to the top of the three arrester units, as shown in

Fig. 6.10 (left).

93

8 (σ̄P ,σ̄P ) 8 (σ̄P ,σ̄I )

235 kV

376 kV

4 470 kV 4 282 kV

423 kV 470 kV

Rel. Error / %

423 kV 329 kV

0 0

0 5 10 15 20 25 30 0 5 10 15 20 25 30

470 kV 329 kV

4 4 282 kV 423 kV

470 kV

423 kV 235 kV 376 kV

0 0

0 5 10 15 20 25 30 0 5 10 15 20 25 30

Nonlinearity Coefficient α

Figure 6.13: Relative error of the total heat dissipation in the resistor column of the

surge arrester, approximated by the sum of the Joulean losses at the

observation points, as function of nonlinearity coefficient α for the four

equivalent conductivity combinations. The different curves correspond

to the peak voltage of the respective simulations. Especially the two

equivalent conductivity combinations on the right-hand side provide

good estimates of total Joulean heat losses for high degrees of material

nonlinearity.

94

its maximum value is not necessarily reached at the topmost observation point.

Therefore, the error of the 470 kV curve for α = 30.89 is about 7 % instead of the

previously cited 3 %

The subsequent Fig. 6.12 shows the absolute value of electric field strength

recorded during a normal, i.e., without envelope equation, electro-quasistatic sim-

ulation of the simplified model of a multi-unit surge arrester for peak voltage

Û = 470 kV and nonlinearity coefficient α = 30.89. Below, the resistive current

density that is calculated from these field strengths is plotted.

The shape of the field strength, and in consequence of the potential, deviates

visibly from a sine shape for the curve obtained at the uppermost observation point

18. There, the peak field strength is reduced due to the varistor material entering

the nonlinear regime. Still, the corresponding points of the lower arrester units

show no significant distortion. For much of the arrester, the assumption of the

envelope equation model, which posits that the potential is essentially harmonic, is

valid, even when the arrester operates close to its threshold voltage. The figure also

shows that there is a minor phase shift between the different observation points.

The reduction of peak field strength for the uppermost observation point is

accompanied by a significant increase of conductivity and resistive current, as

shown for the corresponding curve in the lower part of the figure. The high

nonlinearity of the varistors leads to a sharp rise of the resistive current during a

sub-interval of the voltage period. To put this increase into context, it has to be

noted that the peak capacitive current is still twice as large as the peak resistive

current.

Finally, Fig. 6.13 shows the relative error for the sum of the extracted Joulean heat

losses in the resistor column. The curves prove that the injection of electric energy

into the surge arrester is calculated quite accurately for each of the equivalent

conductivity combinations. However, for high degrees of nonlinearity, the two

combinations on the right-hand side are clearly preferable with relative errors

falling below 1 %.

The second test of the envelope equation model consists of a voltage stress com-

parison. Apart from including a second guard ring, the geometry is the same as

before, but the previously used power-law characteristic is replaced by a charac-

teristic derived from Fig. L.4 of standard IEC 60099-4. The characteristics for σ̄P

and σ̄I are obtained by numerical integration of equations 6.7 and 6.9. The model

corresponds to case C of the standard, a 2D computation with two grading rings and

consideration of resistive effects. The expected maximum voltage of the three units

is given as 41 %/m (top), 39 %/m (middle) and 29 %/m (bottom). With respect to

accuracy, it is stated that [d]eviations of 1 %/m to 2 %/m may typically be expected.

95

45

39.8

40.2

38.9

40

39.3

Voltage Stress / %/m

35

29.2

30

28.8

25

29 ±2 39±2 41 ±2

20

1 2 3

Distance from Arrester Base / m

Figure 6.14: Voltage stress along the resistor column between bottom and top of

the surge arrester. The blue curve is obtained, when a power-equivalent

conductivity is assumed, whereas the violet curve is based on the current-

equivalent conductivity. The numbers in the plot indicate that all results

lie within the range considered as acceptable according to standard

IEC 60099-4 Annex L. However, the brown curve, which represents the

solution of a normal time-domain solution, shows that there is a notable

difference.

96

In Fig. 6.14 the results of the finite element simulations are shown. As heat

generation is of no importance, the figure presents only two solutions for the

envelope-equation model and the solution of a normal transient simulation. The

deviation of envelope-equation model solutions lies well within the the expected

range. It reaches at most 0.8 %/m for the current-equivalent and 1.2 %/m for the

power-equivalent conductivity. Thus, the results of the envelope equation model

satisfy the requirements of standard IEC 60099-4. However, that does not mean

that the envelope equation model provides indeed entirely satisfactory results. The

figure shows that the obtained voltage stress curves differ nonetheless substantially

from the curve obtained by means of a normal transient simulation.

cation of the real problem. This leads inevitably to a loss of accuracy. Given the

difficulty of solving the problem without simplifications, this sacrifice of accuracy

may be tolerated, if it is not too large. To determine the extent of this sacrifice,

many simulations of a surge arrester model have been realized.

While the results of the presented studies confirm the limitations of the model,

the error is not unacceptably high, if the right equivalent material characteristics

are used. Even if the nonlinearity of the material properties and stray capacitances

lead to potentials, which are not perfectly harmonic, the validity of the envelope

equation model for surge arrester simulations is strongly supported by the effects of

field grading and high varistor permittivity, which was set to "r = 800 in accordance

with IEC 60099-4.

The error with respect to heat generation depends strongly on the degree of

nonlinearity α of electrical conductivity. While the results are exact for linear

materials, the error generally reaches its maximum for low values of α and decreases

as the nonlinearity increases further. A possible explanation is that the model leads

to relatively good results, if either the resistive or capacitive currents dominate,

whereas the error is greatest during the transition between these two states.

The total heat generation by Joulean losses in the surge arrester is quite accurately

estimated, if the appropriate equivalent conductivity is chosen. For high values of α

the observed error lies in the range of no more than 1 %. This low level of error can

be seen as a confirmation of the envelope equation model.

The distribution of the heat losses is subject to greater error. Nevertheless, the

maximum observed error of local heat generation does not exceed 10 %. This is

astonishing, given the sensibility of this value with respect to field strength, as heat

generation increases proportionally to Ê α+2 .

97

A second test has been the comparison of the voltage stress along the resistor col-

umn for the envelope equation model with reference values given in a standard and

results of a normal transient field simulation. Given that the degree of nonlinearity

is relatively low, below six, one may expect a larger error following the preceding

observations. Nevertheless, the deviation of the calculated voltage stress from its

reference values lies comfortably inside the range of tolerance. The comparison

with the results of the normal transient simulation shows that the stress curves of

the envelope equation model differ significantly from their expected shape.

Even if the envelope equation model has been applied only to purely electric prob-

lems until now, there are no relevant obstacles to its application for the examination

of the thermal stability of surge arresters. The displayed results are promising,

although it has to be seen to which extent an envelope equation model can be used

for the simulation of surge arresters or other nonlinear problems.

98

7 Nonlinear Circuit Model for

Varistor-Based Stress Control

7.1 Capacitance and Conductance Matrices

The relationship between electric potentials and charges in linear media is described

by a matrix Cinf :

Q = Cinf ϕ. (7.1)

Consequently, the vector of currents I is related to the time derivative of the poten-

tials by:

I = Cinf ϕ̇, (7.2)

and the energy W stored in the electric field is:

W = ϕ T Cinf ϕ. (7.3)

The matrix Cinf that unites these three properties is commonly called capacitance

matrix.

However, the name is ambiguous. In accordance with Lehner [65] and

Kupfmüller [62] the name capacitance matrix will be given to another matrix,

which will be introduced soon. Instead, the name matrix of influence coefficients is

used as name for matrix Cinf . P

Under the assumption ϕ(r) = l ϕl w l (r), the matrix of influence coefficients can

be derived from GAUSS’s law. The charge Q k of conductor k depends on the potential

ϕl of every conductor l (including l = k) [7]:

Z

Qk = D · dA (7.4)

∂Ω

Z

=− " grad ϕ · dA (7.5)

∂Ω

X Z

= ϕl − " grad w l · dA (7.6)

l ∂Ω

X

= ckl ϕl . (7.7)

l

99

The coefficients ckl , which are found in the k-th row and l -th column of the matrix

of influence coefficients, represent some normalized electric flux out of conductor

k due to the potential at conductor l . An important property of this matrix is its

symmetry, as the boundary ∂ Ω can be moved from surrounding conductor k to

surrounding the respective conductor l in the step of Eq. 7.6. This is shown by

Lehner [65] for a homogeneous background material, or more generally by Simonyi

[89].

The matrix of influence coefficients Cinf is closely related to the capacitance matrix

and can be converted easily into capacitance matrix C. While the matrix of influence

coefficients Cinf is characterized by satisfaction of Eq. 7.1, the characteristic property

of capacitance matrix C is the physical meaning of its matrix coefficients. The

individual coefficients can be considered as the capacitance values of capacitors

in an equivalent circuit. While the off-diagonal coefficient Ckl can be seen as the

capacitance of a capacitor connecting the conductors k and l , the entries on the

diagonal Ckk define the self-capacitance of conductor k, i.e., the capacitance between

conductor k and ground.

The charge Q k at conductor k of the equivalent circuit can be calculated by

X

Q k = Ckk ϕk + Ckl (ϕk − ϕl ) . (7.8)

l

(l6=k)

Comparing the previous equation with Eq. 7.1 reveals the relationship between

the coefficients Ckl of the capacitance matrix and the coefficients ckl of the matrix of

influence coefficients and how to convert them:

¨P

m ckm , if k = l

Ckl = , (7.9)

−ckl , if k 6= l

¨P

m Ckm , if k = l

ckl = . (7.10)

−Ckl , if k 6= l

In analogy to Cinf and C, conductance matrices Ginf and G can be introduced. Ma-

trix G will be called conductance matrix and Ginf matrix of conductance coefficients.

The matrix of conductance coefficients Ginf relates a vector of resistive currents I,

which flow out of a conductor, to the potential of the conductors:

I = Ginf ϕ. (7.11)

100

In place of 7.4, one has:

Z

Ik = J · dA (7.12)

∂Ω

Z

=− σ grad ϕ · dA (7.13)

∂Ω

X Z

= ϕl − σ grad w l · dA (7.14)

l ∂Ω

X

= g kl ϕl . (7.15)

l

Similar to the coefficients Ckl , the coefficients of the conductance matrix Gkl

can be interpreted as the conductance values between conductors k and l . The

dependence of the current at conductor k on the coefficients of the conductance

matrix is given by:

X

I k = Gkk ϕk + Gkl (ϕk − ϕl ) . (7.16)

l

(l6=k)

Thus, the conversion rules of Eqs. 7.9 and 7.10 can be used in effectively the same

way to convert between Gkl and g kl .

Thanks to linearity, the matrices can be extracted easily from a few static field

simulations. One possible way is using Eq. 7.1. For a system of n conductors, one

needs at most n simulations of the static electric field with linearly independent

potentials imposed as boundary condition to estimate Cinf . For example, if one sets

the potential of one conductor k to one and the others to zero and calculates the

charges at the conductors by integrating the electric flux density flowing through

their surface, one obtains all values in the k-th column of matrix Cinf , since the value

of the charge is equal to the value in the corresponding row. The matrix Ginf can be

obtained similarly from simulations of the stationary current field.

In cases, when capacitive current Ic and resistive current Ir are not of completely

different magnitude, one has to calculate the total current by using matrix Cinf as

well as Ginf :

I = Ir + Ic (7.17)

= Ginf ϕ + Cinf ϕ̇. (7.18)

This equation allows the calculation of currents I(t) for an externally imposed

evolution of the potential ϕ .

101

7.2 Nonlinear Matrices and Equivalent Circuits

P and the applicability of

the superposition principle, e.g., by assuming Φ(r) = l Φl w l (r).

The superposition principle says that the influence of every source can be de-

termined individually and that the total influence can be obtained by adding the

contributions of the individual sources [92]. If capacitance or conductance is

nonlinear, the assumption is no longer valid.

For example, consider a capacitor with a voltage-dependent capacitance C(U), so

that Q = C(U) U , and another, arbitrary voltage αU . The superposition principle is

in general violated as:

αQ Q

C(αU) 6= = = C(U). (7.19)

αU U

Hence, the coefficients of any matrix relating charges and potentials have to be

functions of the potential of the conductors.

Furthermore, the equations involving currents or energy are no longer fulfilled

at the same time. There is no longer a single matrix that satisfies all the properties

of Eq. 7.1 to 7.3 in the presence of nonlinear materials. The distinction between

Φ1 G11

C11

C21 , G21

Φ2 C31 , G31

C22

Φ3

G33 G22

C32 , G32

C33

and an equivalent circuit representation (right) with variable capacitors

and resistors. The circuit elements represent self- and mutual capaci-

tances, respectively self- and mutual conductances.

102

capacitance and differential capacitance, which has already described in the section

about nonlinear capacitances (Sec. 5.1), has to be taken into account.

While a nonlinear capacitance matrix can not completely substitute its linear

equivalent, one can require that it satisfies at least one of the conditions. In the

following, the nonlinear matrix of influence coefficients is supposed to satisfy the

relationship Q = Cinf ϕ . Similarly, the matrix of conductance coefficients has to

satisfy I = Ginf ϕ .

In analogy to Eq. 7.18, knowledge of these two matrices allows an approxi-

mate calculation of the currents which depends only on the vector of conductor

potentials ϕ :

d

I(t) = Ginf (ϕ(t)) ϕ(t) + Cinf (ϕ(t)) ϕ(t) . (7.20)

dt

As shown, the matrices Cinf (ϕ(t)) and Ginf (ϕ(t)) are convertible into C(ϕ(t)) and

G(ϕ(t)), from which a nonlinear equivalent circuit can be built. Such a nonlinear

equivalent circuit with three conductors is shown in Fig. 7.1.

The following section presents a method to extract the properties of the circuit

elements of an equivalent circuit. The data for the estimation of the properties

are obtained from numerical field simulations, or possibly measurements. After

extraction of the equivalent circuit parameters the matrices of influence coeffi-

cients Cinf (ϕ(t)) and of conductance coefficients Ginf (ϕ(t)) can be assembled and

the currents flowing between the conductors can be calculated.

simulations. This procedure becomes problematic in the presence of nonlinear

materials. The straightforward approach described in Sec. 7.2, which consists of

setting the potential of the individual conductors successively to one or any other

value and observe the resulting charge to assemble the capacitance matrix, fails as

the values of the coefficients depend on potential. While it is possible to repeat the

calculation for different potential values and interpolate the matrix coefficients, this

does not resolve the problem.

The principal reason for the inadequacy of such an approach is the fact that

the capacitive coupling is no direct function of the conductor potentials. The

coupling between two conductors depends principally on the potential difference

and not on the absolute values of the potentials. Therefore, the approach will not

provide satisfactory results for nonlinear problems, unless the potential of all nearby

conductors is close to zero, if one conductor is excited.

103

By contrast, the method, which is described in the following, depends implicitly on

the difference between the conductor potentials. Its starting point is the estimation

of the circuit element parameters of a nonlinear equivalent circuit. In many respects

it is similar to the approach described in Sec. 5.4 for the characterization of materials.

In fact, its principal practical differences are the extension of the previous approach

from scalar to matrix values and its different purpose. Whereas the extraction

of parameters by means of least squares had been presented as an approach to

characterize materials from measurement data, it is now employed to create models

with reduced number of unknowns from field simulations.

The capacitance and conductance of each capacitor or conductor in an equivalent

circuit according to Fig. 7.1 is supposed to depend only on the difference in potential

between their terminal nodes. For example, the capacitance between the conductors

k and l is assumed to be a function of |Φk − Φl |:

required to depend linearly on some unknown parameters Ckli :

X

Ckl (∆Φ) = Ckli w i (∆Φ), (7.22)

i

some basis functions, as in Sec. 5.4. For example, if the capacitance is estimated as

a polynomial of second order, the equation becomes:

0 1 2

Ckl (∆Φ) = Ckl + Ckl ∆Φ + Ckl · ∆Φ2 . (7.23)

This means that the extraction of an entry Ckl of the capacitance matrix requires

the determination of n + 1 unknowns for a polynomial ansatz of order n. As the

equivalent circuit of a multiconductor system with m conductors has m(m + 1)/2

capacitors, one has to determine (n + 1)m(m + 1)/2 unknowns for the capacitance

matrix and an equivalent number of unknowns for the conductance matrix. Thus,

the number of unknowns can easily become large, although it will always be much

lower than for a field simulation.

In the case of the material characterization the parameters are estimated by

minimizing the difference between measured and simulated current densities. In

the present case, the procedure is similar. The simulated total current Isim of the

equivalent circuit at a certain moment t i is the sum of resistive and capacitive

current:

d

Isim (t i ) = Ginf (t i )ϕ(t i ) + Cinf (t)ϕ(t) . (7.24)

dt | t=t i

104

Its value depends on the value of the parameters Ckli and Gkl i

that indirectly

characterize the matrix of influence coefficients Cinf and of conductance coefficients

Ginf , as well as on the vectors ϕ and ϕ̇ that describe the potential and its time

derivative for all conductors. Isim (t i ) is a vector which describes the current flowing

out of each of the conductors at moment t i . This equation can be written in the

form of a system of linear equations with respect to the unknown parameters Ckli

i

and Gkl . By concatenating the systems for a certain number of time moments, e.g.,

p, a matrix with m p rows and (n + 1) m (m + 1)/2 columns is obtained:

Ginf (t 1 )ϕ(t 1 ) + d

(Cinf (t)ϕ(t)) |

Isim (t 1 ) dt t=t 1

.. .. ..

= . (7.25)

. . .

d

Isim (t p ) Ginf (t p )ϕ(t p ) + dt (Cinf (t)ϕ(t)) |

t=t p

However, the parameters are defined as coefficients for the conductors and

capacitors of an equivalent circuit. The current Isim,k flowing out of conductor k at

an arbitrary moment can be moment can be written as:

X

Isim,k =Gkk (|Φk |)Φk + Gkl (|Φk − Φl |) (Φk − Φl )

l

(k6=l)

d X

+ Ckk (|Φk |)Φk + Ckl (|Φk − Φl |) (Φk − Φl ) (7.26)

dt l

(k6=l)

X

=Gkk (|Φk |)Φk + Gkl (|Φk − Φl |) (Φk − Φl )

l

(k6=l)

X

+ Ckk (|Φk |)Φ̇k + Ckl (|Φk − Φl |) (Φ̇k − Φ̇l ) (7.27)

l

(k6=l)

kl k l

+ Φk + (Φk − Φl ). (7.28)

dt l

d t

(k6=l)

The parameters are finally estimated by minimizing the difference between the

simulated total current Isim of the equivalent circuit and the reference solution Iref :

i ,G i

Ckl kl

105

where Iref is a vector of values extracted from numerical field simulations. It assumes

the position that measured current density Jmeas had for the material characterization.

Isim is the product of a matrix A with a vector of the unknown parameters, which

describe the circuit elements:

..

.

k

Gi j

.

minimize
A .. − Iref
. (7.30)

k k

Gi j ,Ci j

Ck

ij

..

.

2

After solution of the least-square problem, the properties of the circuit elements

and the entries of the capacitance matrix C and conductance matrix G can be

calculated for any given set of conductor potentials ϕ . The current flowing between

the conductors is calculated by determining the numerical values of the capacitance

and conductance matrix from present potentials, converting the matrices into matrix

of influence coefficients Cinf and matrix of conductance coefficients Ginf , and finally

using Eq. 7.20.

To apply the described procedure, data from a reference solution Iref are needed,

which were said to be obtained from numerical field simulations. More precisely,

they are obtained from an electro-quasistatic time-domain simulation of the mul-

ticonductor model. The conductors are excited by linearly independent excitation

signals, ideally similar in frequency and amplitude to the signals, for which one

wishes to obtain the equivalent circuit model. At the individual time steps t i , the

potentials ϕl (t i ), the current flowing out of the respective conductor I l (t i ) and

possibly the time-derivative value ϕ̇l (t i ) are recorded to a file. The current I l (t i ) is

calculated by integration over the surface of the conductor l :

Z

∂D

Il = J+ · dA (7.31)

∂ Ωl

∂t

Z

∂

=− (σ(| grad ϕ|) grad ϕ − ("r (| grad ϕ|) "0 grad ϕ) · dA (7.32)

∂Ω

∂t

l

These recorded data contain all the information that is needed to estimate the

parameters of the equivalent circuit model and its associated matrices.

shown for a cable configuration with fictional nonlinear properties. A nonlinear

106

R0 L0

G0 C0

line system, which allows the calculation of voltage and current along the

line by solving the telegrapher’s equation. Application of the equivalent-

circuit characterization approach leads to an estimation of nonlinear

capacitance C 0 and conductance G 0 per unit length.

equivalent-circuit based model with low number of degree of freedom can be derived

that provides results, which are similar to those of the computationally much more

expensive field simulation.

In the case of a cable, the extracted equivalent circuit can be interpreted as

describing the conductance and capacitance per unit length, as represented by

G 0 and C 0 in Fig. 7.2 for a simple two-conductor system. Therefore, the method

provides half of the information needed to simulate the propagation of a current and

Conductor

Conductivity / nS/m

5

Relative Permittivity

4.5 13

4 11

3.5

9

3

2.5 7

PE 0 10 20 30 40

Figure 7.3: Geometry of the studied cable model, a screened PVC control cable and

the fictional material properties of the PVC material. Neither the material

properties nor the excitation is realistic for the given cable, but they serve

to compare the results of a nonlinear equivalent circuit model with those

from a field simulation.

107

Conductor 1 Conductor 4

150 150

100 100

50 50

120 140 160 180 200 120 140 160 180 200

−50 −50

−100 −100

−150 −150

Current per length/ nA/m

Conductor 2 Conductor 5

100

100

50 50

120 140 160 180 200 120 140 160 180 200

−50 −50

−100

−100

Conductor 3

100

Current

50

from FEM simulation

120 140 160 180 200 from equivalent circuit (24 dofs)

−50

−100

Time / ms

simulation and the derived nonlinear equivalent circuit. The capacitors

and conductors of the equivalent circuit are modeled as polynomials of

order three. Under consideration of symmetry, this leads to a model

with only 24 degrees of freedom (dofs). The calculated current using the

capacitance and conductance matrices derived from the nonlinear equiv-

alent circuit largely agrees with the current from the field simulations

and shows the same distorted sine shape. Differences between finite

element and equivalent circuit currents can be observed close to current

peaks.

108

voltage signal along a nonlinear transmission line. However, the other parameters,

R0 and L 0 , have to be obtained otherwise.

For the purpose of demonstrating the method, a model with more than two

wires is used. Its geometry corresponds to the cross section of a common screened

polyvinyl chloride (PVC) control cable, sold under the name (see Fig. 7.3), but the

material properties have been altered.

First, the conductivity of the PVC material has been raised to a level at which

the resistive currents are comparable in magnitude to the capacitive ones at about

50 Hz. Thereby, none of the two matrices is negligible and the estimation of both

matrices makes sense. Second, conductivity and permittivity of the PVC material

have been replaced by second-order polynomials with respect to field strength to

obtain a nonlinear problem (see Fig. 7.3).

The model is simulated by solving the electro-quasistatics equation in time domain.

Each of the conductors is excited by a voltage signal with the same amplitude, but

different initial phase angle and frequency, which differ by some percent from 50 Hz.

The results of this simulation are used to generate a model of reduced order. This

order depends on the number of coupling conductances and capacitances. Since

the simulated cable consists of five wires plus ground, its equivalent circuit has 15

conductors and capacitors. However, this number is reduced to three thanks to

symmetry:

As permittivity and conductivity are quadratic functions of field strength, the ca-

pacitances and conductances are necessarily also at least quadratic. Thus, each

of the three equivalent conductors and capacitors has at least three parameters,

leading to a minimum total number of 18. Since the field between the conductors

is non-homogeneous and the order of the material nonlinearity does not translate

directly to the cable geometry, a higher order might be necessary for the equivalent

conductors and capacitors. In any case, the equivalent circuit model is of much

lower order than the field simulation model.

In Fig. 7.4 the model with nonlinear capacitance and conductance matrix is

compared to the original electro-quasistatic finite element results. The figure shows

the sum of resistive and capacitive current per meter flowing from the individual

conductors to the other conductors and ground. The matrices have been obtained

by evaluation of the finite element results for the interval 0 s to 1 s. However, given

the spread of the conductor excitation frequencies, the current and voltage data

109

Conductor 1

150

Equiv. Circ.

100

Current per length/ nA/m

50 FEM

−50

−100

−150

Time / ms

Figure 7.5: Enlarged view of the current flowing out of the first wire obtained

from an electro-quasistatic field simulation and the derived equivalent

circuit model. Although the two curves are clearly distinct, at least close

to the current peaks, the equivalent circuit model is capable of largely

reproducing the field simulation current.

provide few new information after about 0.2 s or 10 periods. Thus, a much shorter

simulation duration is sufficient for the purpose of extracting the nonlinear matrices.

By contrast, the voltage and current data of the first few milliseconds have to be

excluded, as the initial field distribution influences the quality of the nonlinear

capacitance and conductance matrix in "steady state".

The blue curves represent the current from the finite element simulation. In

combination with the voltage information these curves constitute the input data

for the derivation of an equivalent circuit model with nonlinear capacitance and

conductance matrices. In purple, the currents calculated from these nonlinear

capacitances and conductances are shown. The model consists of six circuit elements,

the three above-mentioned conductors and capacitors, of third polynomial order,

which leads to 24 degrees of freedom. The current calculated with this equivalent

110

circuit model agrees relatively well with the finite element current. Differences can

be seen at the current peaks. Increasing the polynomial degree of the model reduces

the difference, though only by a limited amount. For example, the average of the

Eq

error |I i −I iF E M |/|I iF E M |, with conductor number i , decreases from 5.96 % to 5.85 %,

when the properties of the circuit elements are modeled as fifth-order polynomials

instead of third-order polynomials. Most of the error has to be attributed to the

limitations of the model, as ignoring the influence of third conductors.

In Fig. 7.5 a part of the previous figure is shown enlarged. The purple curve of

the equivalent-circuit based model is quite similar to the blue finite element solution

with some differences close to the current peaks.

111

8 Summary and Outlook

The last chapters of this thesis presented possible solutions to various problems

related to varistors and microvaristors. While the individual sections treat diverse,

independent problems, some of the sections can be considered as forming part of a

chain of actions. The chain reaches from the characterization of a material for the

creation of a numerical model to the generation of a simple equivalent circuit from

numerical field simulations.

The simple, yet apparently unknown least-squares based characterization method

presented in Ch. 5 is a powerful tool to estimate the electrical properties of a

material. Both, conductivity and relative resistivity are estimated simultaneously

as field-strength dependent functions. This characterization method possesses two

great advantages: It requires relatively few measurement data and the shape of

the excitation signal can be (almost) arbitrary. If necessary, an entire characteristic

curve can be extracted directly from a single lightning-current pulse measurement.

The application of this approach to a test sample consisting of microvaristor-

embedded silicone rubber has led to a better understanding of this material. As

variously suspected, the nonlinearity of (micro-)varistor materials includes relative

permittivity. This assumption is strongly supported by the presented results.

Furthermore, the results from a more detailed examination of the data indicate

the presence of some kind of previously unknown memory effect. Section 5.5

describes the first attempts to combine nonlinear material behavior with relaxation.

While the results can hardly be considered satisfactory, they are hopefully of use for

later studies.

Once the material characteristics are known, one can continue with numerical

field simulations. As already mentioned in the introduction of this thesis, the

simulation of a lightning-current pulse injected into a high-voltage surge arrester

has been excluded, although this had initially been the principal objective. By

contrast, the chapter on surge arrester simulations (Ch. 6) contains a short section

on the calculation of a thermally stationary state by iteratively coupling a transient

electro-quasistatic simulation with the static heat conduction problem. By using

this method, a stationary state in thermal equilibrium is obtained. Reassuringly

the same stationary state is reached for different initial conditions. The main

questions concerning this kind of simulation are the importance of knowing the

thermally stationary state, opposed to knowing the limits of thermal stability, and

the accuracy, given the uncertainty about the modeling parameters. Maybe the

112

possibility to calculate the thermally stationary state can be used in conjunction

with measurements to reduce the modeling parameter uncertainty.

The second presented approach for the simulation of surge arresters introduces

for the first time an envelope equation model to electrical power engineering. The

purpose of this innovative approach is to help with the determination of the thermal

stability limits. This important problem is computationally too expensive for a brute

force simulation. While the results are clearly better than expected by the author

and show the potential of this ansatz, they demonstrate also the limitations of this

method. One has to observe, whether the envelope equation can compete with

alternative methods in development.

After the numerical field simulation the results have to be evaluated. Often the

simulations serve to extract capacitance and other matrices, which can then be

inserted into other programs. In Ch. 7 an approach to derive an equivalent circuit

with nonlinear resistors and capacitors from a transient electro-quasistatic field

simulation is presented. This nonlinear equivalent circuit constitutes an intermediate

step towards nonlinear capacitance and conductance matrices. These matrices can

be used as a tool to simulate the propagation of signals on wires with nonlinear

resistive stress control materials or other nonlinear problems.

113

Appendix

Acronyms

ABC absorbing boundary condition

BEM boundary element method

d.c. direct current

DFG Deutsche Forschungsgemeinschaft

dof degree of freedom

EQS electro-quasistatic

FD finite differences (method)

FDTD finite differences in time domain (method)

FEM finite element method

FIT finite integration technique

GRP glass-reinforced plastic

HV high-voltage

IBC impedance boundary condition

IGBT insulated-gate bipolar transistor

MCOV maximum continuous operating voltage

MO metal oxide

MoM method of moments

MOV metal-oxide varistor

NTC negative temperature coefficient

OLS ordinary least squares

PML perfectly matched layer

PVC polyvinyl chloride

RC resistor-capacitor

115

rms root mean squared

SDIRK singly-diagonally implicit Runge-Kutta (method)

SiC silicon carbide

SPD surge-protective device

SVD singular value decomposition

UHV ultra-high voltage

VDR voltage-dependent resistor

ZnO zinc oxide

Lists of Symbols

Notational Conventions

x vector named x

A matrix named A

ẋ time derivative of x

x̃ estimate or prediction of x

x̂ peak value of x

x phasor of x

j imaginary unit

x0 distributed quantity x

x sim quantity x obtained through numerical simulation

x ref quantity x used as reference for comparisons

xT transposed of x

X−1 inverse of X

Electrical Symbols

α nonlinearity coefficient

B magnetic flux density T

C capacitance F

c wave propagation speed m/s

Cd differential capacitance F

Cinf matrix of influence coefficients F

C0 distributed capacitance F/m

C capacitance matrix F

116

D displacement flux density C/m2

E field strength V/m

" permittivity F/m

"ˆ permittivity (per second) of a linear, inert material F/(m s)

"∞ high-frequency limit of dielectric permittivity F/m

"0 electric constant (or vacuum permittivity) F/m

"r relative permittivity 1

"s low-frequency limit of dielectric permittivity F/m

f frequency Hz

G conductance S

Ginf matrix of conductance coefficients (in analogy to Cinf ) S

G0 distributed conductance S/m

G conductance matrix S

H magnetic field strength A/m

I (static) current A

i (time-varying) current A

I vector of currents A

Ic vector of capacitive currents A

Ir vector of resistive currents A

J current density A/m2

Jc capacitive current density A/m2

Jr resistive current density A/m2

λ wave length m

L0 distributed inductance H/m

µ permeability kg m/(A2 s2 )

µ0 permeability of free space H/m

µr relative permeability 1

ω angular frequency 1/s

P polarization density T

pn electric dipole moment (of particle n) Cm

ϕ scalar potential V

Φ phasor potential V

ϕ vector of potentials V

q electric charge C

Q vector of charges C

% charge density C/m3

%0 initial charge density C/m3

R0 distributed resistance Ω/m

σ electrical conductivity S/m

117

σ̄I current-equivalent conductivity (Sec. 6.2) S/m

σ̄P power-equivalent conductivity (Sec. 6.2) S/m

σ̂ conductivity (per second) of a linear, inert material S/(m s)

σ∞ conductivity at high frequency S/m

σs conductivity at low frequency S/m

T period length s

τ characteristic time s

τdiff diffusion time constant s

τ0 relaxation time constant s

τrelax (dielectric) relaxation time constant s

U (static) voltage V

u (time-varying) current A

Uc maximum continuous operating voltage (MCOV) V

Us maximum rms voltage of system V

W energy J

we electric energy density J/m3

wm magnetic energy density J/m3

Thermodynamic Symbols

α diffusivity m2 /s

cp specific heat capacity J/(kg K)

cv volumetric heat capacity J/(m3 K)

ε emissivity 1

h heat transfer coefficient W/m2

λ thermal conductivity W/(K m)

Q̇ (net) heat flux W

q̇ heat flux density W/m2

% mass density kg/m3

σSB Stefan–Boltzmann constant W/(K4 m2 )

ϑ temperature K

ϑ∞ far-field temperature in a gas/fluid K

ϑsur temperature of the surrounding environment K

w̄ time-averaged heat generation W/m3

ẇgen heat generation density W/m3

B matrix associated with boundary conditions

118

K" Stiffness matrix associated with electric permittivity

Kλ Stiffness matrix associated with thermal conductivity

Kσ Stiffness matrix associated with electric conductivity

Mc p % Mass matrix associated with volumetric heat capacity

Ni i -th basis function

θ Factor characterizing the θ -method for time integra-

tion (0 ≤ θ ≤ 1)

r residual

wi i -th weighting function

wV vector associated with volumetric heat generation

A area m2

Bi j i -th B-Spline of order j

d thickness m

δ Dirac operator

∆t length of time interval s

η relative error 1

L characteristic length m

r location vector m

t time s

x vector of unknowns

119

List of Figures

2.1 Photograph of surge arresters. . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2 Cross section of an arrester unit. . . . . . . . . . . . . . . . . . . . . . . . 13

3.2 Some equivalent circuit models for varistors. . . . . . . . . . . . . . . . 31

3.3 Heat capacity and thermal energy of some arrester materials. . . . . . 33

5.2 Photograph of the measurement setup. . . . . . . . . . . . . . . . . . . . 50

5.3 Illustration for the method of chosen points. . . . . . . . . . . . . . . . 52

5.4 B-splines Bik of orders one to three. . . . . . . . . . . . . . . . . . . . . . 56

5.5 The problem with the assumption of constant permittivity. . . . . . . . 59

5.6 Estimated electrical conductivity σ and relative permittivity "r of a

microvaristor sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.7 Measured and predicted current densities. . . . . . . . . . . . . . . . . . 61

5.8 Difference between calculated and measured current density. . . . . . 63

5.9 Convergence of the material characteristic estimates. . . . . . . . . . . 64

5.10 Characteristics for conductivity σ and relative permittivity "r of a

microvaristor sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.11 Measured and predicted current density for the characteristics of Fig.

5.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5.12 Illustration of the error between observed and predicted currents. . . 68

5.13 Electrical conductivity and relative permittivity of a microvaristor-

filled silicone rubber. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.14 Results for the interpolation of material characteristics. . . . . . . . . . 70

5.15 Estimated conductivity σ∞ , relative permittivity "∞ and the charac-

teristic for σs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.2 Flow chart for the calculation of thermally stationary state. . . . . . . 78

6.3 Estimated electrical conductivity of a ZnO varistor. . . . . . . . . . . . 79

6.4 Estimated relative permittivity of a ZnO varistor. . . . . . . . . . . . . . 80

121

6.5 Stationary heat distribution and momentary electric potential of a

surge arrester model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

6.6 Electric field strength along the resistor column of an arrester model. 82

6.7 Temperature along the resistor column of an arrester. . . . . . . . . . . 83

6.8 Extracted conductivity characteristic of a ZnO varistor. . . . . . . . . . 87

6.9 Relative error and Joulean losses of the envelope equation model. . . 88

6.10 Joulean losses at the observation points. . . . . . . . . . . . . . . . . . . 90

6.11 Maximum relative error for each of the four equivalent conductivity

combinations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

6.12 Absolute value of field strength and resistive current density at various

points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

6.13 Relative error of the total heat dissipation in the resistor column. . . 94

6.14 Voltage stress along the resistor column. . . . . . . . . . . . . . . . . . . 96

circuit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

7.2 Equivalent circuit of an elementary component of a two-wire trans-

mission line system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

7.3 Geometry of a cable model and assumed material properties. . . . . . 107

7.4 Comparison of currents from field simulation and equivalent circuit. 108

7.5 Enlarged view of the current flowing out of the first wire. . . . . . . . 110

122

List of Tables

3.1 Maxwell’s equations in differential and integral form. . . . . . . . . . . 17

5.1 Comparison of the attained error level with several material models. 73

123

Bibliography

[1] H. D. Baehr and K. Stephan, Heat and Mass Transfer. Berlin: Springer, 2006.

Wiley, 1989.

conductivity features of ZnO-based varistors using the laser-pulse method,”

Materials Science and Engineering A, vol. 371, no. 1-2, pp. 377–381, Apr.

2004. [Online]. Available: http://www.sciencedirect.com/science/article/

B6TXD-4C09KCG-3/2/b43fbdb98c7c6594107adb9eecba8130

vol. 79, no. 11, pp. 8629–8633, Jun. 1996.

absorption capability of ZnO varistors,” IEEE Transactions on Power Delivery,

vol. 14, no. 1, p. 152–162, Jan. 1999.

Verdrahtungsstrukturen,” Dissertation, Technische Universität Wien, Wien,

Nov. 1994. [Online]. Available: http://www.iue.tuwien.ac.at/phd/bauer/

diss.html

lations based on the finite integration technique and a mixed circuit formula-

tion,” Dissertation, Technische Universität Darmstadt, Darmstadt, 2007.

Society, vol. 82, no. 3, pp. 485–502, 1999. [Online]. Available:

http://dx.doi.org/10.1111/j.1151-2916.1999.tb01793.x

[10] R. W. Clough, “The finite element method in plane stress analysis,” in Proceed-

ings, 2nd Conference on Electronic Computation, A.S.C.E. Structural Division,

Pittsburgh, 1960.

125

[11] R. Courant, “Variational methods for the solution of problems of equilibrium

and vibrations,” Bull. Amer. Math. Soc, vol. 49, no. 1, p. 23, 1943.

McGraw-Hill, 1998.

[13] A. J. Davies, The finite element method : an introduction with partial differen-

tial equations. Oxford; New York: Oxford University Press, 2011.

Aided Geometric Modeling. London: Academic Press, 1993, pp. 27–49.

thermal reaction of zinc oxide to current impulses,” in 30th ICLP 2010,

Cagliari, Italy, 12.-16.9.2010, paper 6B-1331, 2010. [Online]. Available:

http://tubiblio.ulb.tu-darmstadt.de/48583/

simulation of zinc-oxide surge arresters,” in Scientific Computing in Electrical

Engineering SCEE 2010, ser. Mathematics in Industry, B. Michielsen and J.-R.

Poirier, Eds. Springer Berlin Heidelberg, 2012, pp. 213–221. [Online].

Available: http://dx.doi.org/10.1007/978-3-642-22453-9_23

Differentialgleichungen. Berlin: Walter de Gruyter, 2008.

composites for HV applications based on microvaristors,” in Solid Dielectrics,

2004. ICSD 2004. Proceedings of the 2004 IEEE International Conference on,

vol. 1, 2004, pp. 403–406.

grading part 2: Materials and applications,” Electrical Insulation Magazine,

IEEE, vol. 27, no. 2, pp. 18–29, 2011.

[21] K. Eda, “Zinc oxide varistors,” Electrical Insulation Magazine, IEEE, vol. 5,

no. 6, pp. 28–30, 1989. [Online]. Available: 10.1109/57.44606

Journal of Applied Physics, vol. 56, no. 10, pp. 2948–2955, 1984. [Online].

Available: http://link.aip.org/link/?JAP/56/2948/1

126

[23] R. Einzinger, “Metal oxide varistor action – a homojunction breakdown

mechanism,” Applications of Surface Science, vol. 1, no. 3, pp. 329–340, May

1978. [Online]. Available: http://www.sciencedirect.com/science/article/

B6X3T-46RVJJ5-2/2/aae00132aa9fa5169f6c646f77cb97dc

[24] ——, “Metal oxide varistors,” Annual Review of Materials Science, vol. 17,

no. 1, pp. 299–321, 1987. [Online]. Available: http://arjournals.

annualreviews.org/doi/abs/10.1146/annurev.ms.17.080187.001503

Verlag, 1996.

vol. 1.

of numerical analysis, vol. 7, pp. 713–1018, 2000.

in complex materials,” in Fractals, Diffusion, and Relaxation in Disordered

Complex Systems, ser. Advances in Chemical Physics. John Wiley & Sons,

2006, vol. 133.

[29] F. Fernández and R. Diáz, “Metal-oxide surge arrester model for fast transient

simulations,” in Proceedings of the IPST’01, Rio de Janeiro, Jun. 2001.

percolation threshold,” Journal of Applied Polymer Science, vol. 72, no. 12, pp.

1573–1582, 1999. [Online]. Available: http://dx.doi.org/10.1002/(SICI)

1097-4628(19990620)72:12<1573::AID-APP10>3.0.CO;2-6

Fils, 1822.

dans les corps solides,” in Oeuvres de Fourier. Paris: Gauthiers-Villars

et Fils, 1890, vol. 2. [Online]. Available: http://gallica.bnf.fr/ark:

/12148/bpt6k33707/f220

Berlin: Springer, 2007.

127

[34] L. Gaul, M. Kögl, and M. Wagner, Boundary Element Methods for Engineers

and Scientists. An Introductory Course with Advanced Topics. Berlin: Springer,

2003.

tion modeling of sub-cycle dynamics and harmonic generation in nonlinear

waveguides,” Optics express, vol. 15, no. 9, p. 5382–5387, 2007.

[36] G. H. Golub and C. F. Van Loan, Matrix computations, 3rd ed. Baltimore:

Johns Hopkins University Press, 1996.

London; Paris; Tokyo; Hong Kong: Springer, 1990.

Ceramic Society, vol. 73, no. 7, p. 1817–1840, 1990. [Online]. Available:

http://dx.doi.org/10.1111/j.1151-2916.1990.tb05232.x

von MO-Überspannungsableitern durch externe Steuerung,” in Feldsteuernde

Isoliersysteme: Werkstoffe, Design, Prüfung und Simulation ; Vorträge des ETG-

Workshops vom 22. bis 23. November 2011 in Darmstadt. Berlin: VDE-Verlag,

2011.

ser. IEE Power & Energy. London: IEE, 2009, no. 40, pp. 191–244.

Chichester; New York: John Wiley, 1999.

senschaftlichen Rechnens. Wiesbaden: Vieweg + Teubner, 2009.

New Way to Suppress Transients. New York: Electronics, Oct. 1972.

Cliffs: Prentice-Hall (downloadable from Massachusetts Institute of

Technology: MIT OpenCourseWare), 1989, Licence: Creative Commons

Attribution – NonCommercial – Share Alike. [Online]. Available:

http://ocw.mit.edu

128

[45] A. R. Hileman, Insulation Coordination for Power Systems, ser. Power Engi-

neering. New York: Marcel Dekker, 1999, no. 9.

2001.

grading systems for UHV metal-oxide surge arresters a new approach to

numerical simulation and dielectric testing,” in Cigré Session, Paris, Aug.

2008.

5th ed. Boca Raton, Fla.: CRC Press, 2011.

Journal of applied mechanics, vol. 8, no. 4, pp. 169–175, 1941.

Subcommittee Surge Protective Devices, “Modeling of metal oxide surge

arresters,” Power Delivery, IEEE Transactions on, vol. 7, no. 1, pp. 302–309,

Jan. 1992. [Online]. Available: http://ieeexplore.ieee.org/xpls/abs_all.jsp?

arnumber=108922

and Mass Transfer. Singapore: John Wiley & Sons, 2013.

[54] Institution of Electrical Engineers, The lightning flash, ser. IEE power & energy

series. London: Institution of Electrical Engineers, 2003, no. 34.

semiconductor,” Reports on Progress in Physics, vol. 72, no. 12, p. 126501,

2009. [Online]. Available: http://stacks.iop.org/0034-4885/72/i=12/a=

126501

baden: Vieweg + Teubner, 2011.

129

[59] M. Kobayashi, M. Mizuno, T. Aizawa, M. Hayashi, and K. Mitani,

“Development of zinc-oxide non-linear resistors and their applications to

gapless surge arresters,” IEEE Transactions on Power Apparatus and Systems,

vol. PAS-97, no. 4, pp. 1149–1158, Jul. 1978. [Online]. Available: http:

//ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=4181542

in distributed parameter systems by discrete regularization,” Journal of

mathematical analysis and applications, vol. 119, no. 1, p. 128–152, 1986.

3rd ed. Heidelberg: Springer, 2009. [Online]. Available: http:

//site.ebrary.com/id/10313511

Heidelberg: Springer, 2008. [Online]. Available: http://d-nb.info/

989924491/34

Wiesbaden: Vieweg + Teubner, 2011.

[64] M. Lat, “Thermal properties of metal oxide surge arresters,” Power Apparatus

and Systems, IEEE Transactions on, vol. PAS-102, no. 7, pp. 2194–2202, 1983.

Berlin: Springer-Verlag, 2006. [Online]. Available: http://site.ebrary.com/

id/10228862

Verfahren, 4th ed. Berlin, Heidelberg: Springer-Verlag Berlin Heidelberg,

2007.

Journal of Applied Physics, vol. 46, no. 3, p. 1332, 1975. [Online]. Available:

http://link.aip.org/link/JAPIAU/v46/i3/p1332/s1&Agg=doi

Transactions on, vol. 24, no. 2, p. 202–207, 2011.

junction studies of ZnO varistors – theory and experiment,” Applied

Physics Letters, vol. 33, no. 9, p. 830, 1978. [Online]. Available:

http://link.aip.org/link/APPLAB/v33/i9/p830/s1&Agg=doi

130

[70] H.-J. Martin, Die Ferroelektrika, ser. Technisch Physikalische Monographien.

Leipzig: Akademische Verlagsgesellschaft Geest & Portig K.-G., 1964, vol. 15.

of Applied Physics, vol. 10, no. 6, pp. 736–746, Jun. 1971.

nativity and growing process,” Hoden Kenkyu, vol. 48, no. 3, p. 5–17, 2005.

Philosophical Transactions of the Royal Society of London, vol. 155, no. 0, pp.

459–512, Jan. 1865. [Online]. Available: http://rstl.royalsocietypublishing.

org/cgi/doi/10.1098/rstl.1865.0008

[74] D. McHenry, “A lattice analogy for the solution of stress problems,” Journal

of the ICE, vol. 21, no. 2, pp. 59–82, 1943.

E. Hotta, “Multiple beam envelope equations for electron injectors

using a bunch segmentation model,” Physical Review Special Topics -

Accelerators and Beams, vol. 15, no. 6, Jun. 2012. [Online]. Available:

http://link.aps.org/doi/10.1103/PhysRevSTAB.15.064201

[76] M. Modest, Radiative Heat Transfer, ser. Chemical, Petrochemical & Process.

Academic Press, 2003.

tronics, vol. 38, no. 8, pp. 114–117, 1965.

and connections,” Reviews of Geophysics, vol. 37, no. 1, p. 151, 1999. [Online].

Available: http://www.agu.org/pubs/crossref/1999/1998RG900006.shtml

The Computer Journal, pp. 308–313, 1965.

[80] X.-M. Pan, W.-C. Pi, M.-L. Yang, Z. Peng, and X.-Q. Sheng, “Solving problems

with over one billion unknowns by the MLFMA,” Antennas and Propagation,

IEEE Transactions on, vol. 60, no. 5, p. 2571–2574, 2012.

[81] P. Pinceti and M. Giannettoni, “A simplified model for zinc oxide surge

arresters,” Power Delivery, IEEE Transactions on, vol. 14, no. 2, pp. 393–398,

1999. [Online]. Available: 10.1109/61.754079

131

[82] W. Polifke and J. Kopitz, Wärmeübertragung : Grundlagen, analytische und

numerische Methoden. München [u.a.]: Pearson, 2005.

equations,” Applied Numerical Mathematics, Dec. 2012. [Online]. Available:

http://linkinghub.elsevier.com/retrieve/pii/S0168927412002036

nahmevermögen von Metalloxidwiderständen eingesetzt in Hochspan-

nungsnetzen unter Berücksichtigung eines komplexen Fehlerkriteriums,”

Ph.D. dissertation, TU Darmstadt, Aug. 2008.

ergy absorption capability and time-to-failure of varistors used in station-class

metal-oxide surge arresters.” IEEE Transactions on Power Delivery, vol. 12,

no. 1, pp. 203–212, Jan. 1997.

arrester design,” IEEE Transactions on Power Apparatus and Systems,

vol. 96, no. 2, pp. 647–656, Mar. 1977. [Online]. Available: http:

//ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1601978

conditions in electromagnetics. London: Institution of Electrical Engineers,

1995.

waveguide,” IEEE Transactions on Microwave Theory and Techniques,

vol. 17, no. 4, pp. 204–210, Apr. 1969. [Online]. Available: http:

//ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1126932

1st ed. Oxford: Pergamon Press, 1963.

transient electro-quasistatic fields using advanced projection techniques,”

Advances in Radio Science, vol. 4, pp. 49–53, 2006. [Online]. Available:

www.adv-radio-sci.net/4/49/2006/

der Elektrotechnik 1. Wiesbaden: B.G. Teubner Verlag, 2005.

132

[93] R. Süße and B. Marx, Theoretische Elektrotechnik. Band 5: Elektrische Net-

zwerke - Berechnung und Synthese von Schaltungen für vorgegebenes Bifurka-

tionsverhalten. Ilmenau: Wissenschaftsverlag Ilmenau, 2002, vol. 5.

method. Boston: Artech House, 1995.

3rd ed., ser. Artech House antennas and propagation library. Boston: Artech

House, 2005.

junctions” within a ZnO varistor,” Journal of Applied Physics, vol. 61, no. 4, p.

1562, 1987. [Online]. Available: http://link.aip.org/link/JAPIAU/v61/i4/

p1562/s1&Agg=doi

[97] L. T. Tenek and J. H. Argyris, Finite element analysis for composite structures,

ser. Solid mechanics and its applications. Dordrecht ; Boston: Kluwer

Academic Publishers, 1998, no. 59.

for Industrial and Applied Mathematics, 1997.

systems in practical applications. Berlin; New York: Springer, 2001.

voltage field simulations of large scale insulator structures including 2-d

models for nonlinear field-grading material layers,” IEEE Trans. Magn.,

vol. 45, no. 3, pp. 980–983, Mar. 2009. [Online]. Available: http:

//ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=4787325

[101] T. Weiland, “A discretization model for the solution of Maxwell’s equations for

six-component fields,” Archiv Elektronik und Uebertragungstechnik, vol. 31,

pp. 116–120, Mar. 1977.

method to heat conduction analysis,” Nuclear Engineering and Design,

vol. 4, no. 3, pp. 276–286, Oct. 1966. [Online]. Available: http:

//www.sciencedirect.com/science/article/pii/0029549366900513

133

[103] K. S. Yee, “Numerical solution of initial boundary value problems involving

Maxwell’s equations in isotropic media,” IEEE Trans. Antennas and Propaga-

tion, pp. 302–307, 1966.

[104] Q.-H. Zhang and D.-J. Chen, “Percolation threshold and morphology of

composites of conducting carbon black/polypropylene/EVA,” Journal of

Materials Science, vol. 39, no. 5, pp. 1751–1757, Mar. 2004. [Online].

Available: http://dx.doi.org/10.1023/B%3AJMSC.0000016180.42896.0f

thermal stability,” IEEE Trans. Power Delivery, vol. 25, no. 3, pp. 1526–1529,

Jul. 2010. [Online]. Available: http://ieeexplore.ieee.org/lpdocs/epic03/

wrapper.htm?arnumber=5477225

[106] O. C. Zienkewicz, R. Taylor, and J. Zhu, The Finite Element Method. Its Basis

& Fundamentals, 6th ed. Oxford: Elsevier Science & Technology, 2005.

[107] O. C. Zienkiewicz, “The birth of the finite element method and of computa-

tional mechanics,” International Journal for Numerical Methods in Engineering,

vol. 60, no. 1, pp. 3–10, 2004.

problems,” The Engineer, vol. 220, no. 6, pp. 507–510, 1965.

134

Acknowledgements

Lastly, I wish to thank everybody who helped me in one way or another to write this

thesis. My particular thanks go to:

• Prof. Dr.-Ing. Thomas Weiland for his invitation to join TEMF and his role as

thesis supervisor,

• PD Dr. rer. nat. Erion Gjonaj for his immeasurably valuable advice as tutor,

• Max Tuczek, Sébastien Blatt, Jan Debus and Prof. Dr.-Ing. Hinrichsen at the

High Voltage Labs of TU Darmstadt for their contributions and cooperation,

• Prof. Dr.-Ing. Binder and all other colleagues from research group FOR 575,

many others who were willing to give me moral support,

• Dragos Munteanu, Heike Seiler, Achim Veuskens and Marianne Dorn for

ensuring good working conditions,

• all other colleagues for their valuable input, a nice time and many memorable

and great moments.

135