Explicit Time Integration of Transient Eddy Current Problems

TL;DR

This paper introduces an explicit time integration approach for transient eddy current problems by transforming the system into ODEs, utilizing the CFL criterion, and employing PCG with CSPE for efficient computation, validated on a benchmark.

Contribution

It presents a novel explicit integration scheme using a generalized Schur-complement and CSPE to improve efficiency in transient eddy current simulations.

Findings

The scheme is validated successfully on the TEAM 10 benchmark.

Explicit Euler method with CFL stability is feasible for these problems.

The use of CSPE accelerates PCG convergence in the method.

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.

Equations

$$\kappa \frac{\partial \vec{\mathrm{A}}(\mathrm{t})}{\partial \mathrm{t}}+\vec{\nabla }\times \left(\nu \left(\vec{\mathrm{A}}(\mathrm{t})\right)\vec{\nabla }\times \vec{\mathrm{A}}(\mathrm{t})\right)={\vec{\mathrm{J}}}_{\mathrm{s}}(\mathrm{t}),$$

$$\mathbf{M}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{a}+\mathbf{K}(\mathbf{a})\mathbf{a}={\mathbf{j}}_s,$$

$$\left[\begin{array}{cc}{\mathbf{M}}_{cc}& 0\\ 0& 0\end{array}\right]\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{\mathbf{a}}_c\\ {\mathbf{a}}_n\end{array}\right]+\left[\begin{array}{cc}{\mathbf{K}}_{cc}({\mathbf{a}}_c)& {\mathbf{K}}_{cn}\\ {\mathbf{K}}_{cn}^{\mathrm{T}}& {\mathbf{K}}_{nn}\end{array}\right]\left[\begin{array}{c}{\mathbf{a}}_c\\ {\mathbf{a}}_n\end{array}\right]=\left[\begin{array}{c}0\\ {\mathbf{j}}_{s,n}\end{array}\right],$$

$${\mathbf{M}}_{cc}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{a}}_c+\left[{\mathbf{K}}_{cc}({\mathbf{a}}_c)-\left({\mathbf{K}}_{cn}{\mathbf{K}}_{nn}^{\#}{\mathbf{K}}_{cn}^{\mathrm{T}}\right)\right]{\mathbf{a}}_c$$

$${\mathbf{a}}_n$$

$${\mathbf{a}}_c^{m+1}:={\mathbf{a}}_c^m+\Delta t{\mathbf{M}}_{cc}^{-1}\left[{\mathbf{K}}_{cn}{\mathbf{K}}_{nn}^{\#}{\mathbf{j}}_{s,n}^{m+1}-\left({\mathbf{K}}_{cc}({\mathbf{a}}_c^m)-{\mathbf{K}}_{cn}{\mathbf{K}}_{nn}^{\#}{\mathbf{K}}_{cn}^{\mathrm{T}}\right){\mathbf{a}}_c^m\right],$$

$${\mathbf{a}}_n^{m+1}:={\mathbf{K}}_{nn}^{\#}{\mathbf{j}}_{s,n}^{m+1}-{\mathbf{K}}_{nn}^{\#}{\mathbf{K}}_{cn}^{\mathrm{T}}{\mathbf{a}}_c^{m+1}.$$

$${\mathbf{K}}_{nn}{\mathbf{k}}_p$$

$${\mathbf{k}}_p$$

$${\mathbf{V}}^{\mathrm{T}}{\mathbf{K}}_{nn}\mathbf{Vz}={\mathbf{V}}^{\mathrm{T}}\mathbf{r},$$

$${\mathbf{x}}_{0,\mathrm{SPE}}:=\mathbf{Vz},$$

$$\Delta t\le \frac{2}{{\lambda }_{\max }\left({\mathbf{M}}_{cc}^{-1}\left({\mathbf{K}}_{cc}({\mathbf{a}}_c)-{\mathbf{K}}_{cn}{\mathbf{K}}_{nn}^{\#}{\mathbf{K}}_{cn}^{\mathrm{T}}\right)\right)},$$

$${\lambda }_{\max }\left({\mathbf{M}}_{cc}^{-1}\left({\mathbf{K}}_{cc}({\mathbf{a}}_c)-{\mathbf{K}}_{cn}{\mathbf{K}}_{nn}^{\#}{\mathbf{K}}_{cn}^{\mathrm{T}}\right)\right)\propto \frac{1}{h^2\kappa \mu },$$

Full text

\publishers

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.

Explicit Time Integration of Transient Eddy Current Problems

J. Dutiné*∗*

1

  

M. Clemens

1

  

S. Schöps

2

  

G. Wimmer

3

(11affiliationmark: Chair of Electromagnetic Theory, School of Electrical, Information and Media Engineering, University of Wuppertal, Rainer-Gruenter-Str. 21, 42119 Wuppertal, Germany22affiliationmark: Graduate School CE, Technische Universität Darmstadt, Dolivostr. 15, 64293 Darmstadt, Germany33affiliationmark: Fakultät für angewandte Natur- und Geisteswissenschaften, Hochschule für angewandte Wissenschaften Würzburg-Schweinfurt, Ignaz-Schönstr. 11, 97421 Schweinfurt, Germany

) ††footnotetext: 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: 1099-1204). This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Self-Archiving.

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 differential-algebraic 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 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 Newton-Raphson 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 pseudo-inverse 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:

[TABLE]

in which $\kappa$ is the electrical conductivity, $\vec{\mathrm{A}}(\mathrm{t})$ is the time-dependent magnetic vector potential, $\nu$ is the reluctivity which may be nonlinear in ferromagnetic materials and $\vec{\mathrm{J}}_{\mathrm{s}}(\mathrm{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

[TABLE]

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 $\mathbf{j}_{s}$ 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 $\mathbf{a}$ in (2) is reordered into a part $\mathbf{a}_{n}$ corresponding to dofs in the nonconductive regions and into a part $\mathbf{a}_{c}$ corresponding to dofs in conductive regions. The matrices $\mathbf{K}$ and $\mathbf{M}$ are decomposed accordingly and yield the reformulation of (2) with

[TABLE]

where $\mathbf{M}_{cc}$ is the conductivity matrix, $\mathbf{K}_{cc}$ is the part of the curl-curl matrix in conductive media, $\mathbf{K}_{nn}$ is the singular curl-curl part in air and $\mathbf{K}_{cn}$ and is the coupling part of $\mathbf{K}$ between both conducting and non-conducting regions, and $\mathbf{j}_{s,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]

[TABLE]

where $\mathbf{K}^{\#}_{nn}$ is the matrix representation of the pseudo-inverse of $\mathbf{K}_{nn}$. The matrix product $\mathbf{K}_{cn}\mathbf{K}^{\#}_{nn}\mathbf{K}^{\mathrm{T}}_{cn}$ is the generalized Schur complement.

Equation (4a) represents an ODE in which the explicit Euler method with time step width $\Delta t$ can be used. In the $(m+1)$-th time step

[TABLE]

is computed. With $\mathbf{a}^{m+1}_{c}$ the solution vector $\mathbf{a}^{m+1}_{n}$ for the degrees of freedom in the nonconductive region can be separately calculated with

[TABLE]

The singular matrix $\mathbf{K}_{nn}$ 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:

[TABLE]

where $\mathbf{r}$ represents one of the vectors with which $\mathbf{K}^{\#}_{nn}$ is multiplied in (5), (6), i.e. $\mathbf{j}^{m+1}_{s,n}$, $\mathbf{K}^{\mathrm{T}}_{cn}\mathbf{a}^{m}_{c}$ and $\mathbf{K}^{\mathrm{T}}_{cn}\mathbf{a}^{m+1}_{c}$. The matrix $\mathbf{K}^{\#}_{nn}$ does not need to be computed explicitely. Alternatively, a vector $\mathbf{k}_{p}$ according to (7a), (7b) is computed by the PCG method for each right-hand side $\mathbf{r}$ stated above replacing a matrix vector multiplication $\mathbf{K}^{\#}_{nn}\mathbf{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 $\mathbf{M}_{cc}$ in (5) and (6) algebraic systems of equations are solved by the PCG method. Due to constant matrices $\mathbf{M}_{cc}$ and $\mathbf{K}_{nn}$ the corresponding systems of equations form multiple right-hand side (mrhs) problems.

Solutions for $\mathbf{k}_{p}$ 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 $\mathbf{V}=\{\mathbf{v}_{1}\mid...\mid\mathbf{v}_{m}\}$. Solving the projected system:

[TABLE]

where $\mathbf{r}$ is the varying right-hand side vector, with a direct method yields the vector $\mathbf{z}\in\mathbb{R}^{m}$. This holds the coefficients for computing a new start vector

[TABLE]

by linear combination of the earlier computed column vectors of the operator $\mathbf{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 $\mathbf{V}$ by MGS, and is referred to as $\mathbf{v^{m+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, $\mathbf{v^{m+1}}$ is appended as column vector to the operator $\mathbf{V}$. Otherwise, it is inserted into the last column of the operator $\mathbf{V}$. An increasing number of PCG iterations results from increasingly worse start vectors obtained with an operator $\mathbf{V}$ storing insufficient information with respect to the new vector $\mathbf{r}$. Appending $\mathbf{v^{m+1}}$ increases the spectral information content of $\mathbf{V}$ without deleting any information from previously included column vectors. In the computation of the product $\mathbf{K}_{nn}\mathbf{V}$ in (8) only the matrix-column-vector product $\mathbf{K}_{nn}$ $\mathbf{v^{m+1}}$ changes in every time step. The other matrix-vector products $\mathbf{K}_{nn}$ $\mathbf{v_{i}},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 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/\tau))$ 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

[TABLE]

as stated in [1]. For the maximum eigenvalue $\lambda_{\mathrm{max}}$ the proportionality:

[TABLE]

holds, whereas $\mu$ is the permeability, $\kappa$ is the conductivity and $h$ is the smallest edge length of the mesh used for spatial discretization.

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 simulation time. The maximum stable time step width for this mesh is $\Delta t=1.2\,\mu s$. As stated in Section 3, the number of column vectors in the CSPE operator $\mathbf{V}$ is increased if the number of PCG iterations required in each time step is increasing. Preventing the number of column vectors in $\mathbf{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 ($N_{CG,acc.}$). The number of column vectors in V is only increased if $N_{CG,acc.}$ is exceeded. Parameters varied in the simulations are $N_{CG,acc.}$ and the tolerance for the stopping criterion for the PCG method when computing the pseudo-inverse. For each PCG tolerance $tol_{PCG}\in\{10^{-8},10^{-7},10^{-6}\}$ simulations with $N_{CG,acc.}\in\{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 3 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 $1^{\mathrm{st}}$ order edge elements.

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 5 and 5. For each PCG tolerance, the use of the CSPE method reduces the average number of PCG iterations compared to using the solution 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 3 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.

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-curl-matrix 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.