Survey on Semi-Explicit Time Integration of Eddy Current Problems

09/20/2017 ∙ by Jennifer Dutiné, et al. ∙ 0

The spatial discretization of the magnetic vector potential formulation of magnetoquasistatic field problems results in an infinitely stiff differential-algebraic equation system. It is transformed into a finitely stiff ordinary differential equation system by applying a generalized Schur complement. Applying the explicit Euler time integration scheme to this system results in a small maximum stable time step size. Fast computations are required in every time step to yield an acceptable overall simulation time. Several acceleration methods are presented.



There are no comments yet.


page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Spatially discretizing the magnetic vector potential formulation of eddy current problems, e.g by the Finite Element Method (FEM), yields a differential-algebraic equation system (DAE) [1]. Commonly, only unconditionally stable implicit time integration methods as e.g. the implicit Euler method or the singly diagonal implicit Runge-Kutta schemes can be used for time integration of the infinitely stiff DAE system [2]. In every time step at least one large nonlinear algebraic equation system needs to be solved due to the nonlinear BH-characteristic in ferromagnetic materials. The Newton-Raphson method is frequently used for linearization and requires at least one iteration per time step. Here, the Jacobian and the resulting stiffness matrix are updated in each iteration and the resulting linear algebraic equation system needs to be solved efficiently.

Applying explicit time integration schemes avoids the necessity of linearization, because nonlinearities only appear in right-hand side expressions. A first approach to use an explicit time integration method for eddy current problems has been proposed in [3], where in the conducting regions of the problems the Finite Difference Time Domain (FDTD) method is used. In the nonconducting regions, i.e., in the air, the corresponding parts of the solution are computed using the boundary element method (BEM) [3]. In a second approach presented in [4], the Discontinuous Galerkin FEM and an explicit time integration method are used for computations in the conducting regions. Continuous FEM ansatz functions and an implicit time integration scheme are applied to the nonconducting regions of the problem [4]. Both approaches in [3] and [4] are based on a separate treatment of conducting and nonconducting regions. A different approach presented in [1] and [5] proposes a Schur complement reformulation of the eddy current problem. In[6] the use of a generalized Schur complement is proposed. Here, a pseudo-inverse of the singular curl-curl matrix in nonconducting regions is evaluated using the preconditioned conjugate gradient (PCG) method. This evaluation forms a multiple-right hand side problem and suitable start vectors for the PCG method are computed using the cascaded Subspace Projection Extrapolation (CSPE) method, which is a modification of the Subspace Projection Extrapolation (SPE) method [7, 6]. Alternatively, the Proper Orthogonal Decomposition (POD) method can be used for computing improved start vectors [8]. Computations can be accelerated further by using a selective update strategy for updating the reluctivity matrix in conducting regions [9]. This paper presents a survey on the methods presented in [6, 8, 9].

2 Mathematical Formulation

The partial differential equation


describes magnetoquasistatic field problems using the time-dependent magnetic vector potential , where is the electrical conductivity, is the eventually ferromagnetic, i.e., nonlinearly field dependent, reluctivity and is the transient source current density.

Discretizing (1) in space, e.g. by FEM using edge elements [10, 11], yields a differential-algebraic equation system (DAE) described by


where is the mass-matrix, is the time dependent vector of the magnetic vector potential, is the stiffness-matrix and

is the vector of the transient source currents. The degrees of freedom (DoFs) are separated into two vectors

and for conducting and nonconducting material, respectively and (2) is re-ordered into


where is the conductivity matrix in conducting regions, is the nonlinear part of the reluctivity related stiffness matrix in conducting regions, is the part of the curl-curl operator in air, which is singular, and is the source current vector corresponding to excitations in nonconducting regions. is positive definite if using a conventional Galerkin scheme with (possibly high-order) edge elements as test and ansatz functions [10, 11].

The generalized Schur complement expression


where is the matrix representation of a pseudo-inverse of , is applied to (3) and transforms the DAE into an ordinary differential equation (ODE) system


for the vector , i.e., the degrees of freedom only situated in conductive material [5, 1, 6]. The preconditioned conjugate gradient (PCG) method is used for evaluating a pseudo-inverse of [6]. Alternatively, the singular matrix can be regularized using a grad-div regularization by which is transformed into the discrete Laplacian operator in free space [5]. Due to finite stiffness, (5) can be integrated in time using explicit time integration schemes as e.g. the explicit Euler method. Here, in the -th time step the expressions


are evaluated for the degrees of freedom in the conductor domain and in the nonconductive domains consecutively, where is the time step size.

Evaluating a pseudo-inverse of and the inverse of in (7) and (8) repeatedly using the PCG method forms multiple right-hand side (mrhs) problems since the matrices involved remain constant. The subspace extrapolation (SPE) method can be used for computing improved start vectors for the PCG method [7, 6]. Solution vectors from previous time steps are orthonormalized using the modified Gram-Schmidt method and form the linearly independent column vectors of the operator . The projected system


where represents the new right-hand side for the full system, is solved for using a direct method. The linear combination of the column vectors in weighted with the coefficients in yields the improved start vector :


Only the last column vector in the operator changes in every time step. Therefore, when computing in (9), all other matrix-column-vector products evaluated can be reused from previous time steps. This modification of the SPE start vector generation method is referred to as ”cascaded SPE” (CSPE).

Alternatively, the proper orthogonal decomposition (POD) method can be used for computing improved start vectors for the PCG method [8, 12]. A snapshot matrix

is assembled using solutions from previous time steps as column vectors. This matrix is decomposed by the singular value decomposition (SVD)

[13] into:


where is a diagonal matrix of the singular values and and are orthogonal matrices. The first column vectors of corresponding to the largest singular values , for which holds


become the column vectors of the reduced matrix with a threshold value . A threshold value commonly used in practical computations is . The improved start vector for the PCG method is computed by


The explicit Euler method is only stable for time step sizes smaller than a Courant-Friedrich-Levy-type time step size given by [1]:


where the maximum eigenvalue

is proportional to


assuming non-singularity of . Here, is the smallest edge length in the mesh, is the electrical conductivity, and is the permeability. Numerical tests have shown that

unfortunately does not give a sharp estimate of

, such that the largest eigenvalue has to be computed numerically.

Fine spatial discretization can result in small stable time step sizes, due to (15), that can be in the micro- to nano second range. Considering the dynamics of the usual excitation currents in magnetoquasistatic problems, this corresponds to a strong over-sampling. It is assumed that the excitation current does not change significantly between succeeding time steps. Therefore it is expected that the vector in (7), (8) also only changes marginally between succeeding time steps. The matrix is thus only updated if the change between the vector at the time step and the vector from the time step at which the matrix has last been updated is larger than a chosen tolerance , as described by [9]:


where denotes the l2-norm. However, depending on the gauging used, a different norm might be more appropriate, e.g. using the magnetic energy norm.

(a) TEAM 10 model geometry and position S1. Steel plates are colored in blue and red, the coil in green. There is a 0.5 mm wide air gap between the blue and red steel plates.
(b) Comparison of simulation results using a mesh of 700,000 DoFs and the measured results published in [14] at position S1.
Figure 1: TEAM 10 model geometry and comparison of simulation results.

3 Numerical Validation

The ferromagnetic TEAM 10 benchmark problem is spatially discretized using first order edge element FEM ansatz functions [14, 15]. The model geometry is shown in Fig. 0(a). The excitation current is described by a function. A time interval of duration is calculated. The accuracy of the employed simulation code is proven using an implicit time integration method and a fine mesh discretization of about 700,000 DoFs. The resulting average magnetic flux density is compared with the measurement results published in [14] in Fig. 0(b). As this simulation takes a simulation time of 5.38 days on a workstation with an Intel Xeon E5 processor, a coarser mesh is applied for the simulations using the explicit Euler method for time integration. The applied coarse spatial discretization yields 29,532 DoFs and results in a maximum stable time step size , such that 100,000 explicit Euler time steps are required for this problem.

Computing improved start vectors for the PCG method using either CSPE or POD reduces the average number of required PCG iterations compared to using the solution from the previous time step as start vector. An algebraic multigrid method is used as preconditioner. The results for the evaluation of the pseudo-inverse of using a PCG tolerance of are shown in Fig. 1(a). Using the selective update strategy for updating the matrix does not significantly decrease accuracy, as is shown in Fig. 1(d). The number of required updates and the simulation time are significantly reduced, as is depicted in Fig. 1(c) and Fig. 1(b). If is updated in every time step 100,000 updates are performed during the entire simulation. A workstation with an Intel Xeon E5 processor and an NVIDIA TESLA K80 GPU are used for these simulations. The matrix is inverted directly using GPU acceleration. This is only possible, as the matrix is only of dimension 5955x5955 in this test problem. For more refined discretizations the PCG method should be used for inverting the matrix .

(a) Averagely required number of PCG iterations for evaluating the pseudo-inverse of using either CSPE, POD, or the solution from the previous time step as start vector for the PCG method.
(b) Comparison of simulation times for a simulation using implicit Euler method, the explicit Euler method with updates of in every time step and the explicit Euler method using the selective update strategy and different tolerances .
(c) Number of updates of for different tolerances in (17).
(d) Average magnetic flux density at position S1.
Figure 2: Numerical results for implicit and various variants of the explicit Euler method.

4 Conclusion

The application of a generalized Schur complement to the spatially discretized magnetic vector potential formulation of magnetoquasistatic field problems transformed a DAE of infinite stiffness into a finitely stiff system of ODEs. This ODE system is integrated with the explicit Euler method. For the evaluation of a pseudo-inverse the PCG method was adopted. Improved start vectors were computed with the CSPE and the POD method, reducing the number of required PCG iterations in simulations of the ferromagnetic TEAM 10 benchmark problem. A selective update strategy for the reluctivity matrix taking into account the specific problem dynamics reduced the number of required updates and the simulation time. So far, the small stable time step size of the explicit Euler method results in high computational effort which can be overcome using massive GPU-parallelization to reduce the required computational time per time step significantly.


This work was supported by the German Research Foundation DFG (grant numbers CL143/11-1 and SCHO1562/1-1). The third author is supported by the Excellence Initiative of the German Federal and State Governments and The Graduate School of Computational Engineering at TU Darmstadt.


  • [1] Schöps, S., Bartel, A., Clemens, M.: Higher order half-explicit time integration of eddy current problems using domain substructuring. IEEE Trans. Mag., 48(2), 623–626 (2012)
  • [2] Hairer, E., Wanner, G.: Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd rev. edn.). Springer, Berlin (1996)
  • [3] Yioultsis, T.V., Charitou, K.S., Antonopoulos, C.S., Tsiboukis, T.D.: A finite difference time domain scheme for transient eddy current problems. IEEE Trans. Mag., 37(5), 3145–3149 (2001)
  • [4] Außerhofer, S., Bíro, O., Preis, K.: Discontinuous galerkin finite elements in time domain eddy-current problems. IEEE Trans. Mag., 45(3), 1300–1303 (2009)
  • [5] Clemens, M., Schöps, S., De Gersem, H., Bartel, A.: Decomposition and regularization of nonlinear anisotropic curl-curl daes. Compel, 30(6), 1701–1714 (2011)
  • [6] Dutiné, J., Clemens, M., Schöps, S., Wimmer, G.: Explicit time integration of transient eddy current problems. presented at the 10th International Symposium on Electric and Magnetic Fields (EMF2016), full paper submitted.
  • [7] Clemens, M., Wilke, M., Schuhmann, R., Weiland, T.: Subspace projection extrapolation scheme for transient field simulations. IEEE Trans. Mag., 40(2), 934–937 (2004)
  • [8] Dutiné, J., Clemens, M., Schöps, S.: Multiple right-hand side techniques in semi-explicit time integration methods for transient eddy current problems, accepted for presentation at the 17th Biennial IEEE Conference on Electromagnetic Field Computation (CEFC 2016).
  • [9] Dutiné, J., Clemens, M., Schöps, S.: Explicit Time Integration of Eddy Current Problems Using a Selective Matrix Update Strategy. accepted for presentation at the 17th International IGTE Symposium on Numerical Field Calculation in Electrical Engineering (IGTE 2016).
  • [10] Monk, P.: Finite element methods for Maxwell’s equations. Oxford University Press, Oxford (2003)
  • [11] Jin, J.-M.: The finite element method in electromagnetics (3rd Edition), Wiley-IEEE Press, Hoboken (2014)
  • [12] Chatterjee, A.: An introduction to the proper orthogonal decomposition. Current Sci., 78(7), 808–817 (2000)
  • [13] Trefethen, L.N., Bau, D.: Numerical Linear Algebra. Society for Industrial and Applied Mathematics, Philadelphia (1997)
  • [14] Nakata, T., Fujiwara, K.: Results for benchmark problem 10 (steel plates around a coil). Compel, 9(3), 181–192 (1990)
  • [15] Kameari, A.: Calculation of transient 3d eddy current using edge-elements. IEEE Trans. Mag., 26(2), 466–469 (1990)