A time-filtering method for plasma simulation: High bulk conductivity

doi:10.1016/j.apm.2021.02.017

Applied Mathematical Modelling

Volume 95, July 2021, Pages 447-467

https://doi.org/10.1016/j.apm.2021.02.017 Get rights and content

Highlights

•
Stable and efficient solver for dielectric barrier discharges with multiple materials.
•
The time step is not constrained by large variations of charge mobilities in different solid layers.
•
Applications to realistic electrical components affected by internal partial discharges.

Abstract

In this paper we aim to make possible the simulation of plasma - solid interaction in real insulating components using a set of first-principle partial differential equations. The high ratio between the conductivities of the many materials present in electrical components causes the associated numerical problem to become very stiff and this curtails the maximum allowed time step increase.

To simulate some realistic components using sufficiently large time steps a filtering method is applied. A novel theoretical analysis has been developed showing that the filter allows the use of larger time steps without affecting the accuracy of the method. This analysis is backed by the results of a significant numerical experiment.

Introduction

In this work, we aim to develop a proper filtering time-integration scheme for the simulation of the interaction between plasma and a dielectric bulk where the conductivity of some of its sub-regions is relatively high.

The simulation capabilities of a breakdown plasma at ambient pressure have quickly evolved in the last years: a first example is the simulation of the corona phenomenon using a set of drift-diffusion equations coupled with an electrostatic equation [1], [2]. The model then evolved to include the interaction of plasma with dielectric surfaces and is also known as dielectric barrier discharges (DBDs), see [3], [4].

There are many technical applications linked to this field however, in this work, we underline that these numerical tools have been very useful to simulate the evolution of the aging of electrical components. In particular, in many cases the aging is linked to the partial discharges (PDs) phenomenon, see [5]: a sequence of discharges takes place in internal voids inside a dielectric bulk causing a progressive de-polymerization of the material leading to a degeneration of its insulating performances. This process is also known as treeing [6], [7].

The simulation, based on a set of partial differential equations (PDEs), has proven to be a valuable tool in this particular branch of electrical engineering even though the considered geometries are, by now, limited to some very simple cases [8], [9], [10] that are far more simplified than those of real electrical components [11]. In contrast to previous literature, here we want to extend the capabilities developed by authors in [1], [10] to include some real configurations. However, a major difficulty has arisen: in many electrical components, some semi-conductive materials are used to avoid any local increase of the electric field that, in turn, may trigger the formation of an electrical tree. These materials can be characterized by very high conductivity and this has a huge impact on the simulation. In spite of the fact the algorithms used in [1], [10] are fully implicit and asymptotic preserving [12], [13], [14], the resulting nonlinear system of equations, that derives from the discretization of the continuous model, is so stiff that it cannot be solved except by constraining the maximum time step increase to unacceptable values. This happens because the time scales associated with the movement of charges in some sub-regions of the bulk are very small with respect to other time scales present in the model.

The proper numerical handling of problems with different time scales has been treated in various fields of scientific computing. For instance, the simulation of atmospheric phenomena [15], [16] provides an example: fluids dynamics in domains that exhibit vertical scales quite different from the horizontal ones is characterized very different time scales [17], [18], [19]. A similar problem is observed in the simulation of compressible fluids dynamics. At low speeds, the acoustic waves travel much faster than the transport of mass thus generating two very different time scales, see, for instance, [20], [21]. Moreover, the inclusion of complex chemical models into atmospheric dynamic simulation codes poses another difficulty. The evolution of chemical species is ruled, in many cases, by sets of complex, highly nonlinear, differential equations that embed a large number of time scales that cannot be easily decoupled from each other. Similar problems have also been observed in the simulation of combustion processes [22], [23]. Other relevant fields are the simulation of shallow water equations [24], [25] and the simulation of turbulence where the time scales involved in the evolution of tiny turbulent structures are very different from the ones associated with the mean drift of the fluid. Finally, even the simulation of plasma discharges in a gas provides a case where several time scales are present. For instance, plasma relaxation time during discharge events can be very small and the model must be treated with some specific numerical schemes such as the asymptotic preserving ones [12], [13], [26].

In almost all the cases described above some specific solutions have been developed that involve a proper manipulation of the original set of equations or a proper time discretization. Many general-purpose algorithms can also be found in the literature provided that the system can be cast as a canonical set of ordinary differential equations (ODEs). This is a rather common case since, using some proper space discretization schemes, a PDE, such as the one discussed in [27], [28], [29], [30], can be reduced to an ODE. In [31], [32], [33], [34], [35], [36] a rather large review of the methods for stiff ODEs can be found. Besides these general-purpose approaches, some specific ones have been developed. For instance, [37] addresses the problem of a rapid initial transitory phase which enforces the use of small time steps. In that work, this problem is circumvented by properly preparing the initial data to avoid the reconstruction of the transitory phase. In [38] the case of a fast oscillatory problem coupled with a slow one is analyzed. It is also treated how to avoid a too severe time step restriction while conserving the accuracy of the reconstruction of the slow component. Also time-splitting algorithms [39], [40] are quite popular and effective especially in all cases where the slow and fast components can be easily decoupled. Some high order schemes for stiff equations, based on the Taylor series expansion, have been discussed in [41]. Finally, a benchmark of the performances of various stiff solvers can be found in [42].

Almost all the methods proposed for stiff ODEs use some implicit discretization techniques and this implies that a set of non-linear equations must be solved for each time step. The resulting set of non-linear equations is, in many cases, as stiff as it was the original problem and, even though many specific solvers [43], [44] have been developed, their solution is equally difficult. This means that non-linear solvers may fail to converge triggering a step rejection and the reduction of the time increment. Fortunately, there is a class of time integration schemes where the original ODE is filtered so that the corresponding nonlinear problem is easier to solve or, in other words, is characterized by a lower Lipschitz constant. There are many filtering methods and they have been developed using different starting points. In [45] it is pointed out that even an explicit scheme can be used effectively using large time steps without incurring instabilities. The key point is to use two different time increments. Steps using large time increments are interposed with sets of steps characterized by small time increments. If large and small time steps are chosen effectively, the number of small time steps required to obtain a stable scheme depends only mildly on the Lipschitz constant, i.e. on its logarithm. A very similar idea is exploited in [46], [47] where a multi-scale approach is used. A set of steps with a small time increment is used to estimate the source term of the steps using larger time increments. A slightly different approach is employed in [48], [49] where a dynamic filter is coupled to the original ODE to eliminate the fast and undesired time scales.

Multi-step approaches [45], [46] and [47] require a very precise knowledge of the time scales involved, which is a rather difficult task when dealing with the discretization of complex sets of PDE. Moreover, since plasma models usually include the concentration of chemical species, we must ensure that all concentrations remain positive during the simulation runs. In general, it is not straightforward to guarantee the positivity of the solution using the latter two methods. On the contrary, method described in [48], [49] seems to be adaptable to our case and we have chosen it for our applications.

The main contribution of this work is to adapt the [48], [49] technique to make possible the simulation of internal discharges in electrical components using a PDE-based, first principles, set of equations detailed in [10], [50]. Moreover, we will also show that the proposed technique preserves the accuracy while greatly reducing the computational cost.

Let us now briefly review the structure of the paper: in Section 2 we introduce a complete plasma-bulk interaction model that embeds both the plasma dynamics described in [1], [2], [10] and the movement of charges in the bulk as described in [4]. Since this complete plasma model is, from a mathematical point of view, nearly un-tractable we will present, in Section 3, a reduced model embedding the presence of a very high bulk conductivity coefficient. We will also introduce a proper time semi-discretization and a fixed-point scheme for the solution of the resulting non-linear problem as treated in [51], [52]. The resulting iterative solution of a sequence of linear problems is tackled with a proper finite-element and finite-volume space discretization. In Section 4 we analyze the iterative scheme and prove that, if the bulk conductivity coefficient is large, its performances are poor. A proper fix is introduced in Section 5 and its properties are analyzed in Section 6 from a theoretical point of view. The implementation of the fix for the original complete plasma-bulk model is briefly discussed in Section 7 and its results are verified by a technically relevant numerical experiment in Section 8. A critical review of the results obtained is included in Section 9.

Section snippets

Plasma model

Let $Ω$ be an open domain partitioned into two non-overlapping regions: $Ω^{g}$ and $Ω^{s}$ such that $\bar{Ω} = {\bar{Ω}}^{g} \cup {\bar{Ω}}^{s}$ . The domain $Ω^{g}$ is occupied by a gas (usually air) while the domain $Ω^{s}$ is occupied by an insulating medium. This volume can be partitioned into $N^{s}$ non-overlapping subregions: $Ω^{s, α}$ with $α = 1, \dots, N^{s}$ and ${\bar{Ω}}^{s} = \cup_{α = 1}^{N^{s}} {\bar{Ω}}^{s, α}$ . We denote with $\partial Ω^{g},$ $\partial Ω^{s, α}$ the boundaries of the domains $Ω^{g}$ and $Ω^{s, α}$ respectively and let $Γ^{g} = \partial Ω^{s} \cap \partial Ω^{g},$ $Γ^{g, α} = \partial Ω^{s, α} \cap \partial Ω^{g}$ and $Γ^{s, α, β} = \partial Ω^{s, α} \cap \partial Ω^{s, β}$ with $α, β = 1, \dots, N^{s}$ be the surfaces separating

Reduced bulk conduction model and its discretization

The problem described by Eqs. (1), (2) and (3) is nearly intractable therefore we introduce here a reduced problem along with a simplified nomenclature. The original problem will be considered once again in the results section. Let $Ω$ be an open domain, then we consider the following dimensionless problem describing the drift of charges in a solid medium ${\begin{matrix} \frac{\partial u}{\partial t} + ν \vec{\nabla} \cdot (u \vec{E}) = 0, & in Ω, \\ \vec{E} = - \vec{\nabla} ϕ, & in Ω, \\ \vec{\nabla} \cdot \vec{E} = \frac{λ}{ν} u, & in Ω, \\ ϕ = ϕ_{b}, & on \partial Ω, \\ u = u_{b}, & on \partial Ω, \end{matrix}$ to be solved for $t > 0$ where a proper set of initial conditions $u = u^{0}$ has been

Fixed point convergence analysis

In this section we analyze the convergence properties of the iterative scheme (11) using the approach already introduced in [12], [13], [14], [51] to analyze the asymptotic properties of implicit and semi-implicit time schemes. This will involve a slight modification of the proposed problem to aid the analytical treatment: for instance an infinite domain is considered. We set $u^{(p)} = 1 + δ u^{(p)},$ ${\vec{E}}^{(p)} = \hat{E} + δ {\vec{E}}^{(p)}$ where $\hat{E}$ is a constant unit vector $| \hat{E} | = 1$ and $δ u^{(p)} ≪ 1,$ $| δ {\vec{E}}^{(p)} | ≪ 1$ o.e. in $Ω$ . Substituting

Filtering method

We have outlined how the iterative scheme (11) may show a slow convergence if the bulk mobility coefficient is high. To avoid this problem we apply the strategy outlined in [48], [49] complementing the system with a filter. In this case we consider a simple, first order dynamic system: $\frac{d u}{d t} = - \frac{1}{τ_{a}} (u - v),$ where $u$ is the output, $v$ is the input signal and $τ_{a}$ is a characteristic cut-off time. We discretize (29) to obtain the following implicit scheme: $u^{n + 1} - u^{n} = - \frac{Δ t}{τ_{a}} u^{n + 1} + \frac{Δ t}{τ_{a}} v^{n + 1} .$ Now, rearranging the terms

Accuracy analysis

In this section we study how the introduction of the $θ$ parameter modifies the solution of (11) and, in particular, we will consider fully discrete schemes. We start by studying the difference between the discrete solution associated with the scheme described in (11) ${\begin{matrix} \vec{\nabla} \cdot ((1 + λ^{t} u^{(p)}) \vec{E}) = \frac{λ^{t}}{ν^{t}} u^{n}, & in Ω, \\ \vec{E} = - \vec{\nabla} ϕ, & in Ω, \\ u - u^{n} + ν^{t} \vec{\nabla} \cdot (u \vec{E}) = 0, & in Ω, \\ ϕ = ϕ_{b}, & on \partial Ω, \\ u = u_{b}, & on \partial Ω, \end{matrix}$ where we have dropped, for the sake of clarity, the iteration index $(p + 1),$ and the modified system ${\begin{matrix} \vec{\nabla} \cdot ((1 + θ λ^{t} u_{θ}^{(p)}) {\vec{E}}_{θ}) = \frac{λ^{t}}{ν^{t}} u^{n}, & in Ω, \\ {\vec{E}}_{θ} = - \vec{\nabla} ϕ_{θ}, & in Ω, \\ u_{θ} - u^{n} \end{matrix}$

Implementation details

The technique presented in Section 5 and analyzed in Section 6 has been applied to the model described in Section 2. In particular, the problem lies in Eqs. (1) where coefficients $μ_{e}^{s, α},$ $μ_{p}^{s, α}$ may be large for some $α \in [1, N^{s}]$ . Therefore, we introduce a coefficient $θ^{s, α} \in R$ for each $Ω^{s, α},$ $α = 1, \dots, N^{s}$ and we substitute (1) with the following system characterized by a reduced mobility ${\begin{matrix} \frac{\partial n_{e}^{s, α}}{\partial t} - \vec{\nabla} \cdot (n_{e}^{s, α} θ^{s, α} μ_{e}^{s, α} {\vec{E}}^{s, α}) = C_{e}^{s, α}, \\ \frac{\partial n_{p}^{s, α}}{\partial t} + \vec{\nabla} \cdot (n_{p}^{s, α} θ^{s, α} μ_{p}^{s, α} {\vec{E}}^{s, α}) = C_{p}^{s, α}, \\ \frac{d n_{e t}^{s, α}}{d t} = C_{e t}^{s, α}, \\ \frac{d n_{p t}^{s, α}}{d t} = C_{p t}^{s, α}, \\ \frac{\partial}{\partial t} (ρ^{s, α} \end{matrix}$

Numeric experiment

We consider here a technically relevant case that has been the main driver for the development of the methods discussed in this paper. We consider a realistic high voltage cable joint with an internal gaseous defect. Five different solid subdomains are present therefore $N^{s} = 5$ . The first domain $Ω^{s, 1},$ outlined in red in Fig. 2, represents the insulating part of the joint and is usually made of flexible silicon materials. Domains $Ω^{s, 2}$ and $Ω^{s, 3}$ are the insulating parts of the two cables and are, in

Conclusions

In this work we have shown that the filtering technique described in [48], [49], when applied to a charge drift problem, is equivalent solving the same problem but using some reduced mobility coefficients. Although it is not intuitive, we have shown that there are conditions under which the solution of the modified problem does not differ too much from the original one. In particular, this happens when the $λ_{h}^{t}$ and $ν_{h}^{t}$ coefficients are high. The higher they are, the less the solution of the

Acknowledgment

This work has been financed by the Research Found for the Italian Electrical System under the Contract Agreement between RSE and the Ministry of Economic Development. The authors wish to thank L. Barbareschi for her valuable contribution and suggestions.

References (66)

A. Villa et al.
An efficient algorithm for corona simulation with complex chemical models
J. Comput. Phys.
(2017)
A. Villa et al.
A PDE-based partial discharge simulator
J. Comput. Phys.
(2017)
F. Cordier et al.
An asymptotic-preserving all-speed scheme for the euler and navier stokes equations
J. Comput. Phys.
(2012)
P. Degond et al.
Asymptotic-preserving particle-in-cell methods for the vlasov-maxwell system near quasi-neutrality
J. Comput. Phys.
(2017)
A. Sandu et al.
Benchmarking stiff ODE solvers for atmospheric chemistry problems i: implicit versus explicit
Atmos. Environ.
(1997)
A. Sandu
Benchmarking stiff ODE solvers for atmospheric chemistry problems II: rosenbrock solvers
Atmos. Environ.
(1997)
P. Bartello
A comparison of time discretization schemes for two-timescale problems in geophysical fluid dynamics
J. Comput. Phys.
(2002)
F. Cordier et al.
An asymptotic-preserving all-speed scheme for the euler and navier stokes equations
J. Comput. Phys.
(2012)
A. Imren et al.
On the merits of extrapolation-based stiff ODE solvers for combustion CFD
Combust. Flame
(2016)
T. Dubois et al.
Multilevel schemes for the shallow water equations
J. Comput. Phys.
(2005)

R.M.M. Mattheij

On computing smooth solutions of problems with large lipschitz constants

Appl. Numer. Math.

(1986)

G.D. Byrne

Stiff ODE solvers: a review of corning attractions

J. Comput. Phys.

(1987)

B. Sportisse

An analysis of operator splitting techniques in the stiff case

J. Comput. Phys.

(2000)

A. Villa et al.

An symptotic preserving scheme fo the streamer simulation

J. Comput. Phys.

(2013)

A. Villa et al.

An implicit three-dimensional fractional step method for the simulation of the corona phenomenon

Appl. Math. Comput.

(2017)

P. Causin et al.

A dual-mixed hybrid formulation for fluid mechanics problems: mathematical analysis and application to semiconductor process technology

Comput. Methods Appl. Mech. Eng.

(2003)

A. Villa et al.

Three dimensional simulation of the dynamics of electro active polymers using shell elements

Applied Mathematics Computation

(2020)

S. Dolce et al.

Test instrument for the automatic compliance check of cast resin insulated windings for power transformers

Measurement

(2017)

T.N. Tran et al.

Numerical modelling of negative discharges in air with experimental validation

J. Phys. D Appl. Phys.

(2011)

W.S. Kang et al.

Numerical study on the influence of barrier arrangements of dielectric barrier discharge characteristics

IEEE T. Plasma Sci.

(2003)

Y. Serdyuk et al.

Computer modeling of interaction of gas discharge plasma with solid dielectric barriers

IEEE Trans. Dielectr. Electr. Insul.

(2005)

I.W. McAllister

Partial discharges in spheroidal voids. void orientation

IEEE Trans. Dielectr. Electr. Insul.

(1997)

S. Rowland et al.

Application of FEA to image-based models of electrical trees with uniform conductivity

IEEE Trans. Dielectr. Electr. Insul.

(2015)

J.M. Rodriguez-Serna et al.

Numerical simulation of temperature and pressure changes due to partial discharges in spherical cavities within solid dielectrics at different ageing conditions

Energies

(2019)

G. Callender et al.

Critical analysis of partial discharge dynamics in air filled spherical voids

J. Phys. D Appl. Phys.

(2018)

G. Callender et al.

Simulating partial discharge activity in a cylindrical void using a model of plasma dynamics

J. Phys. D Appl. Phys.

(2019)

D. Fabiani et al.

HVDC Cable design and space charge accumulation. part 3: effect of temperature gradient

IEEE Electr. Insul. Mag.

(2008)

J. Schtz

An asymptotic preserving method for linear systems of balance laws based on galerkins method

J. Sci. Comput.

(2014)

G.L. Browning et al.

Analysis of periodic updating for systems with multiple timescales

J. Atmos. Sci.

(1996)

A.G. Peddle

Components of nonlinear oscillation and optimal averaging for stiff PDEs

(2018)

S. Noelle et al.

A weakly asymptotic preserving low mach number scheme for the euler equations of gas dynamics

Computational Methods in Science and Engineering

(2012)

S. Weiqi et al.

A multi-timescale and correlated dynamic adaptive chemistry and transport (CO-DACT) method for computationally efficient modeling of jet fuel combustion with detailed chemistry and transport

Combust. Flame

(2017)

B.T. Nadiga et al.

On simulating flows with multiple time scales using a method of averages

Theor. Comput. Fluid Dyn.

(1997)

Cited by (2)

A multi-resolution method for internal partial discharge simulation
2023, Journal of Computational Physics
Partial discharges and the electrical treeing represent one of key aspects for the ageing of the insulating parts of electrical components. To predict the end of life of dielectric systems several different modelling approaches have been developed. One of the most promising is the modelling using first principles and a set of partial differential equations. In this field, one of the toughest challenges is to cope with the very different scales of the treeing (each branch is a few micrometers wide) and of electrical components (raging from a few centimetres to a few meters).
In this paper we propose an approach, based on the finite element and finite volume methods, which uses two completely independent grids to describe internal defects and electrical components. Some interpolation algorithms have been developed to couple the solution on the two grids. We have also tested the accuracy of our approach by comparing it with previously published methods that use only one grid.
Simulation of surface-plasma interaction with high surface conductivity
2022, Journal of Computational Physics
Citation Excerpt :
This difference between the two solutions also contributes to the fact that the second run has required more steps to reach the same final simulation time. In this paper we have analyzed the problem associated with high values of surface mobility and we have shown that, unfortunately, this is a tougher problem than the bulk conductivity problem even when some high bulk mobility coefficients are considered, see [18]. We have developed an algorithm to deal with this problem and we have shown, both from a theoretically and experimental point of view, that this approach can diminish the total computational cost while preserving the accuracy.
Plasma simulation is getting increasingly important to reproduce technically relevant configurations in electrical engineering. For instance, simulation tools are used to represent the evolution of partial discharges in internal defects such as voids and electrical treeing. In these cases the simulation includes both the gaseous defect and the polymeric bulk where the defect is embedded. A major role is played by the physical characteristics of the surface separating the defect from the bulk and, among these parameters, the surface conductivity is one of the most important aspects. The latter may take very high values due to the cumulative interaction of dielectric materials with a large number of partial discharges.
The simulation of configurations involving high surface conductivities is particularly tough since, as we will show, this may reduce the maximum time step increase of simulation algorithms. In this work, we introduce a proper numerical filter that can circumvent this problem and we show that, using this technique, we can significantly enhance the performances of simulation algorithms without compromising their accuracy. The performances of this approach will be shown in a set of numerical examples.

View full text

A time-filtering method for plasma simulation: High bulk conductivity

Highlights

Abstract

Introduction

Section snippets

Plasma model

Reduced bulk conduction model and its discretization

Fixed point convergence analysis

Filtering method

Accuracy analysis

Implementation details

Numeric experiment

Conclusions

Acknowledgment

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

Atmos. Environ.

Atmos. Environ.

J. Comput. Phys.

J. Comput. Phys.

Combust. Flame

J. Comput. Phys.

Appl. Numer. Math.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

Appl. Math. Comput.

Comput. Methods Appl. Mech. Eng.

Applied Mathematics Computation

Measurement

Numerical modelling of negative discharges in air with experimental validation

J. Phys. D Appl. Phys.

Numerical study on the influence of barrier arrangements of dielectric barrier discharge characteristics

IEEE T. Plasma Sci.

Computer modeling of interaction of gas discharge plasma with solid dielectric barriers

IEEE Trans. Dielectr. Electr. Insul.

Partial discharges in spheroidal voids. void orientation

IEEE Trans. Dielectr. Electr. Insul.

Application of FEA to image-based models of electrical trees with uniform conductivity

IEEE Trans. Dielectr. Electr. Insul.

Numerical simulation of temperature and pressure changes due to partial discharges in spherical cavities within solid dielectrics at different ageing conditions

Energies

Critical analysis of partial discharge dynamics in air filled spherical voids

J. Phys. D Appl. Phys.

Simulating partial discharge activity in a cylindrical void using a model of plasma dynamics

J. Phys. D Appl. Phys.

HVDC Cable design and space charge accumulation. part 3: effect of temperature gradient

IEEE Electr. Insul. Mag.

An asymptotic preserving method for linear systems of balance laws based on galerkins method

J. Sci. Comput.

Analysis of periodic updating for systems with multiple timescales

J. Atmos. Sci.

Components of nonlinear oscillation and optimal averaging for stiff PDEs

A weakly asymptotic preserving low mach number scheme for the euler equations of gas dynamics

Computational Methods in Science and Engineering

A multi-timescale and correlated dynamic adaptive chemistry and transport (CO-DACT) method for computationally efficient modeling of jet fuel combustion with detailed chemistry and transport

Combust. Flame

On simulating flows with multiple time scales using a method of averages

Theor. Comput. Fluid Dyn.