Explicit Time Integration of Transient Eddy
Current Problems

arXiv:1701.03009v1 [cs.CE] 11 Jan 2017

J. Dutiné1∗ , M. Clemens1 , S. Schöps2 , G. Wimmer3
1 Chair of Electromagnetic Theory, School of Electrical, Information and Media Engineering,
University of Wuppertal, Rainer-Gruenter-Str. 21, 42119 Wuppertal, Germany
2 Graduate School CE, Technische Universität Darmstadt, Dolivostr. 15, 64293 Darmstadt,
Germany
3 Fakultät für angewandte Natur- und Geisteswissenschaften, Hochschule für angewandte
Wissenschaften Würzburg-Schweinfurt, Ignaz-Schönstr. 11, 97421 Schweinfurt, Germany
Abstract. For time integration of transient eddy current problems commonly implicit time integration methods are used, where in every time step
one or several nonlinear systems of equations have to be linearized with the
Newton-Raphson method due to ferromagnetic materials involved. In this
paper, a generalized Schur-complement is applied to the magnetic vector potential formulation, which converts a differential-algebraic equation system of
index 1 into a system of ordinary differential equations (ODE) with reduced
stiffness. For the time integration of this ODE system of equations, the explicit Euler method is applied. The Courant-Friedrich-Levy (CFL) stability
criterion of explicit time integration methods may result in small time steps.
Applying a pseudo-inverse of the discrete curl-curl operator in nonconducting
regions of the problem is required in every time step. For the computation of
the pseudo-inverse, the preconditioned conjugate gradient (PCG) method is
used. The cascaded Subspace Extrapolation method (CSPE) is presented to
produce suitable start vectors for these PCG iterations. The resulting scheme
is validated using the TEAM 10 benchmark problem.

1 Introduction
In the computation of magnetoquasistatic fields spatial discretization e.g by the Finite
Element Method (FEM) of the magnetic vector potential formulation yields a differentialalgebraic equation system of index 1 (DAE(1)) [1]. This DAE(1) system is of infinite
stiffness and can only be integrated in time by suitable, unconditionally stable implicit
This is the pre-peer reviewed version of the article, which has been accepted for publication by the
International Journal of Numerical Modelling: Electronic Networks, Devices and Fields (ISSN: 10991204). This article may be used for non-commercial purposes in accordance with Wiley Terms and
Conditions for Self-Archiving.

1

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

time integration methods. Frequently used implicit time integration methods are e.g.
the implicit Euler method or the singly diagonal implicit Runge-Kutta schemes (SDIRK)
[2]. The nonlinear BH-characteristic in ferromagnetic materials requires the solution of
at least one large nonlinear equation system in every time step. For this, the NewtonRaphson method is an established linearization scheme which may require several iterations per time step. In each Newton-Raphson iteration the stiffness-matrix and the
Jacobian-matrix need to be reassembled or at least updated and the resulting linear algebraic equation system has to be solved. As these systems may be high dimensional for 3D
problems, each time step is rather time consuming and the possibilities for acceleration
and parallelization involve complicated numerical linear algebra techniques.
The use of the Newton-Raphson method can be avoided within explicit time integration schemes which are not standard schemes for eddy current problems. An explicit
time integration approach has been first proposed in [3], where the explicit Finite Difference Time Domain (FDTD) method is used in the conductive regions of the problem.
Here, in the nonconductive regions the boundary element method (BEM) is applied
for computing the air parts of the solution, corresponding to a magnetostatic Poisson
problem [3]. As an alternative to the BEM, the use of an adapted Perfectly Matched
Layer (PML) is proposed in [4]. In [5] and [6], the conductive and nonconductive regions
are also treated differently. Here, within a Discontinuous Galerkin FEM-framework the
explicit time integration method is used for regions with conductive materials, while
an FEM formulation with continuous ansatz functions and an implicit time integration
method are applied to the nonconductive regions [5],[6].
The work presented in this paper is originally based on an approach proposed in [1] and
[7] based on a Schur-complement reformulation of the magnetoquasistatic problem. In
this paper, the use of a generalized Schur-complement is proposed in which a pseudoinverse of the singular curl-curl matrix is considered in the nonconductive regions. The
solution of the resulting system of equations is obtained by the preconditioned conjugate gradient (PCG) method. As this results in a multiple right-hand side problem, an
optimized starting vector for the PCG-method can be computed by the subspace projection extrapolation method (SPE) [8]. The algorithm of the SPE has been modified
to a cascaded SPE (CSPE) scheme that is used throughout this work. The generalized
Schur-complement formulation will be explained in Section 2, SPE and CSPE will be
outlined in Section 3 and the corresponding numerical results will be presented in Section
4.

2 Mathematical Formulation
The partial differential equation describing maqnetoquasistatic problems using the magnetic vector potential is given by:
κ

 


~
∂ A(t)
~ × ν A(t)
~
~ × A(t)
~
+∇
∇
= ~Js (t),
∂t

2

(1)

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

~
in which κ is the electrical conductivity, A(t)
is the time-dependent magnetic vector
potential, ν is the reluctivity which may be nonlinear in ferromagnetic materials and
~Js (t) is the time-dependent source current density.
Spatial discretization of (1) by e.g. FEM [9] or the finite integration technique (FIT)
[10], results in a differential-algebraic system of equations of index 1 (DAE(1)) of the
form
d
a + K(a)a = js ,
(2)
dt
where M is the mass-matrix of conductivities, a is the vector of the degrees of freedom
(dofs) corresponding to the magnetic vector potential, K is the stiffness-matrix of reluctivities corresponding to the singular curl-curl operator discretization and js is the
transient source current density. A solenoidal right-hand side in (2) ensures the existence
of a solution, which is achieved if the source current density equals the curl of the electric
vector potential [11].
The dofs vector a in (2) is reordered into a part an corresponding to dofs in the
nonconductive regions and into a part ac corresponding to dofs in conductive regions.
The matrices K and M are decomposed accordingly and yield the reformulation of (2)
with
   

   
Mcc 0 d ac
Kcc (ac ) Kcn ac
0
+
=
,
(3)
0
0 dt an
KT
K
a
j
nn
n
s,n
cn
M

where Mcc is the conductivity matrix, Kcc is the part of the curl-curl matrix in conductive media, Knn is the singular curl-curl part in air and Kcn and is the coupling part of
K between both conducting and non-conducting regions, and js,n is the source current
density in the nonconductive region. Applying the Schur-complement to (3), transforms
the infinitely stiff DAE(1) into an ordinary differential equation system (ODE) of finite
stiffness [7]

Mcc

h

i
d
T
ac + Kcc (ac ) − Kcn K#
K
ac = −Kcn K−1
nn cn
nn js,n ,
dt
#
T
an = K#
nn js,n − Knn Kcn ac ,

(4a)
(4b)

where K#
nn is the matrix representation of the pseudo-inverse of Knn . The matrix
T
product Kcn K#
nn Kcn is the generalized Schur complement.
Equation (4a) represents an ODE in which the explicit Euler method with time step
width ∆t can be used. In the (m + 1)-th time step
h

 i
m
−1
# m+1
m
#
T
m
:=
am+1
a
+
∆tM
K
K
j
−
K
(a
)
−
K
K
K
(5)
cn nn s,n
cc c
cn nn cn ac ,
c
c
cc
is computed. With am+1
the solution vector am+1
for the degrees of freedom in the
c
n
nonconductive region can be separately calculated with
m+1
#
T m+1
:= K#
am+1
.
n
nn js,n − Knn Kcn ac

3

(6)

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

The singular matrix Knn could be regularized with a grad-div-regularization resulting
in a discrete vector Laplacian operator in free space [7]. As this is computationally
expensive [7], here the computation of a pseudo-inverse is proposed using the PCG
method. In (5) and (6), the pseudo-inverse is evaluated by solving systems of the kind:
Knn kp
kp

= r,
:= K#
nn r,

(7a)
(7b)

where r represents one of the vectors with which K#
nn is multiplied in (5), (6), i.e.
#
m+1
T
m
T
m+1
js,n , Kcn ac and Kcn ac . The matrix Knn does not need to be computed explicitely.
Alternatively, a vector kp according to (7a), (7b) is computed by the PCG method for
each right-hand side r stated above replacing a matrix vector multiplication K#
nn r in
(5), (6) with the discrete pseudo-inverse operator.

3 SPE and CSPE
For the evaluation of the pseudo-inverse and for the application of the inverse of the
matrix Mcc in (5) and (6) algebraic systems of equations are solved by the PCG method.
Due to constant matrices Mcc and Knn the corresponding systems of equations form
multiple right-hand side (mrhs) problems.
Solutions for kp from previous time steps are used to obtain an improved start vector
for the PCG method. Within the subspace projection extrapolation (SPE) start vector
generation method [8], the solution vectors from m previous time steps are orthonormalized by the modified Gram-Schmidt process (MGS) to form the linearly independent
column vectors of an operator V = {v1 | ... | vm }. Solving the projected system:
VT Knn Vz = VT r,

(8)

where r is the varying right-hand side vector, with a direct method yields the vector
z ∈ Rm . This holds the coefficients for computing a new start vector
x0,SPE := Vz,

(9)

by linear combination of the earlier computed column vectors of the operator V [8].
In this work, the SPE method is employed as follows: In the (m + 1)-th time step,
the solution from the last, i.e. m-th, time step is orthonormalized against the column
vectors in V by MGS, and is referred to as vm+1 . If the number of required PCG
iterations is larger in the m-th time step than it has been in the (m − 1)-th time step,
vm+1 is appended as column vector to the operator V. Otherwise, it is inserted into
the last column of the operator V. An increasing number of PCG iterations results
from increasingly worse start vectors obtained with an operator V storing insufficient
information with respect to the new vector r. Appending vm+1 increases the spectral
information content of V without deleting any information from previously included
column vectors. In the computation of the product Knn V in (8) only the matrixcolumn-vector product Knn vm+1 changes in every time step. The other matrix-vector
products Knn vi , i = 1, ..., m, have been computed at previous time steps and can be
reused. This modified version of the SPE is the so called ”Cascaded SPE” (CSPE).

4

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

4 Numerical Results
The nonlinear magnetoquasistatic TEAM benchmark problem ”TEAM 10” [12] is simulated. It consists of two steel plates in the shape of square brackets that are placed
facing each other with a third, rectangular steel plate in their midst, with small air gaps
in between, as depicted in Figure 1. At time t = 0 s, the current in the excitation coil
starts to increase according to a (1 − exp(−t/τ )) function and the reaction of the field
of the magnetic flux density was simulated for the first 120 ms. The time step width for
the explicit Euler method is bounded by
2

∆t ≤
λmax



M−1
cc



T
Kcc (ac ) − Kcn K#
nn Kcn

 ,

(10)

as stated in [1]. For the maximum eigenvalue λmax the proportionality:



T
λmax M−1
Kcc (ac ) − Kcn K#
∝
cc
nn Kcn

1
,
h2 κµ

(11)

holds, whereas µ is the permeability, κ is the conductivity and h is the smallest edge
length of the mesh used for spatial discretization.

Figure 1: Geometry of the TEAM 10 model. The line S1 qualitatively shows where the
average magnetic flux density is measured
Eqn. (11) shows that a fine mesh resolution, e.g. of an narrow air gap, results in a
small maximum stable CFL-time step criterion. A rather coarse mesh with 29,532 dofs is
employed in order to allow for the variation of several parameters in an acceptable overall

5

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

simulation time. The maximum stable time step width for this mesh is ∆t = 1.2 µs. As
stated in Section 3, the number of column vectors in the CSPE operator V is increased
if the number of PCG iterations required in each time step is increasing. Preventing
the number of column vectors in V from becoming too large, thus resulting in more
expensive computations in the evaluation of (8), is done by additionally setting a limit
of accepted PCG iterations (NCG,acc. ). The number of column vectors in V is only
increased if NCG,acc. is exceeded. Parameters varied in the simulations are NCG,acc.
and the tolerance for the stopping criterion for the PCG method when computing the
pseudo-inverse. For each PCG tolerance tolP CG ∈ {10−8 , 10−7 , 10−6 } simulations with
NCG,acc. ∈ {1, 3, 5} are run, in order to check the effect on the accuracy of the calculated
magnetic flux density and on the simulation time.
Such a coarse mesh does not yield sufficient accuray with respect to the measured
reference results published in [12]. Therefore, the results for the magnetic flux density
obtained in simulations with the explicit time integration scheme are compared with the
results of a reference simulation on the same mesh using the standard formulation and
the implicit Euler method for time integration. The results are depcited in Figure 2
and show good agreement. In order to prove the accuracy of the employed code itself
using an implicit time integration method, a finer discretization with about 700,000 dofs
is chosen, resulting in a simulation time of 5.38 days on a workstation with an Intel
Xeon E5-2660 processor. The results are presented in Figure 3. Independent of the time
integration method, the spatial disretization is done by FEM based on 1st order edge
elements.

Figure 2: Comparison of results of implicit Figure 3: Comparison of simulation with
and explicit time integration at
about 700,000 dofs and measured
position S1
results stated in [12]
The number of averagely used PCG iterations is decreased by reducing the prescribed
PCG tolerance and by increasing the number of column vectors in the CSPE operator
V, as becomes evident in Figures 4 and 5. For each PCG tolerance, the use of the CSPE
method reduces the average number of PCG iterations compared to using the solution

6

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

Figure 4: Average number of required PCG Figure 5: Maximum number of column veciterations
tors in the CSPE operator V

Figure 6: Duration of the simulations with explicit time integration
from the previous time step as start vector by a factor of at least 4 up to a maximum
of a factor 12. The effect of the number of column vectors in the operator V on the
average number of required PCG iterations is less pronounced for lower PCG tolerances.
Figure 2 demonstrates that there is no loss in accuracy by using a PCG tolerance of 10−6
when computing the pseudo-inverse for this test example. Figure 5 shows that a rather
small number of column vectors of less than 20 in the CSPE operater V is sufficient.
The simulation time is reduced by relaxing the PCG tolerance and using CSPE, as is
shown in Figure 6. A decrease of simulation time with these two schemes by a factor of
about 2.22 was achieved. However, the time needed for simulating this problem with an
implicit Euler method was 2.35 h, so another speed-up of at least a factor of 1.9 will be
necessary for the approach proposed in this method to become faster than the standard
method using an implicit time integration scheme.

7

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

5 Conclusion
A generalized Schur-complement was applied to the spatially discretized magnetic vector
potential formulaiton of the eddy current problem. This transformed an inifitely stiff
DAE(1) into an ODE of finite stiffness that was integrated in time using the explicit
Euler method, thus avoiding linearization by the Newton-Raphson method. Within
the generalized Schur-complement approach a pseudo-inverse for the singular curl-curlmatrix was evaluated by the PCG method. The average number of PCG iterations for
evaluating this pseudo-inverse could be significantly reduced by generating an improved
start vector for PCG with the CSPE method, which was demonstrated on simulations
of the ferromagnetic TEAM 10 benchmark problem. Although the simulation time was
reduced by use of the CSPE method, the main objective of being faster than the implicit
time integration method should only become visible in case of easier to accomplish
parallel implementation for massive many core systems. Further investigations should
focus on further increasing the speed with which each time step is executed and on
reducing the number of required time steps by using coarse space discretization combined
with higher order schemes.

Acknowledgement
This work was supported by DFG grants CL-143/11-1 and SCHO-1562/1- 1, the Excellence Initiative of German Federal and State Governments and the Graduate School CE
at TU Darmstadt.

References
[1] Schöps S, Bartel A, Clemens M. Higher Order Half-Explicit Time Integration of
Eddy Current Problems Using Domain Substructuring. IEEE Transactions on Magnetics 2012; 48(2), pp. 623-626, DOI:10.1109/TMAG.2011.2172780.
[2] Hairer E, Wanner G. 2010. Solving Ordinary Differential Equations II: Stiff and
Differential-Algebraic Problems (2nd rev. edn). Springer-Verlag.
[3] Yioultsis TV, Charitou KS, Antonopoulos CS, Tsiboukis TD. A Finite Difference
Time Domain Scheme for Transient Eddy Current Problems. IEEE Transactions on
Magnetics 2001; 37(5), pp. 3145-3149, DOI:10.1109/20.952563.
[4] Yioultsis TV. Explicit Finite-Difference Time-Domain Method With Perfectly
Matched Layers for Transient Low-Frequency Eddy-Current Analysis. IEEE Transactions on Magnetics 2002; 38(5), pp. 3439-3447, DOI:10.1109/TMAG.2002.802742.
[5] Außerhofer S, Bı́ró O, Preis K. Discontinuous Galerkin Formulation for Eddy-Current
Problems. COMPEL 2009; 28(4), pp. 1081-1090, DOI:10.1108/03321640910959125.

8

J. Dutiné, et al. – Explicit Time Integration of Transient Eddy Current Problems

[6] Außerhofer S, Bı́ró O, Preis K. Discontinuous Galerkin Finite Elements in Time
Domain Eddy-Current Problems. IEEE Transactions on Magnetics 2009; 45(3),
pp. 1300-1303, DOI:10.1109/TMAG.2009.2012604.
[7] Clemens M, Schöps S, De Gersem H, Bartel A. Decomposition and Regularization
of Nonlinear Anisotropic Curl-Curl DAEs. COMPEL 2011; 30(6), pp. 1701-1714,
DOI:10.1108/03321641111168039.
[8] Clemens M, Wilke M, Schuhmann R, Weiland T. Subspace Projection Extrapolation Scheme for Transient Field Simulations. IEEE Transactions on Magnetics 2004; 40(2), pp. 934-937, DOI:10.1109/TMAG.2004.824583.
[9] Kameari A. Calculation of Transient 3D Eddy Current Using Edge-Elements. IEEE
Transactions on Magnetics 1990; 26(5), pp. 466-469, DOI:10.1109/20.106354.
[10] Clemens M, Weiland T. Transient Eddy Current Caclulation with the
FI-Method. IEEE Transactions on Magnetics 1999; 35(3), pp. 1163-1166,
DOI:10.1109/20.767155.
[11] Ren Z. Influence of the R.H.S. on the Convergence Behaviour of the CurlCurl Equation. IEEE Transactions on Magnetics 1996; 32(3), pp. 655-658,
DOI:10.1109/20.497323.
[12] Nakata T, Fujiwara K. Results for Benchmark Problem 10 (Steel Plates Around a
Coil). COMPEL 1990; 9(3), pp. 181-192, DOI:10.1108/eb010074.

9

