Semiexplicit time integration of a reduced magnetic vector potential magneto-quasistatic field formulation

In this article, a reduced vector potential formulation for transient magneto-quasistatic field problems is presented. The finite element method for spatial discretization yields a system of differential algebraic equations (DAE). Applying a generalized Schur complement the DAE system decomposes into a system of ordinary differential equations (ODE) and an algebraic equation. In order to avoid computationally expensive linearization schemes to solve nonlinear eddy current problems the ODE system is integrated in time using a semiexplicit Euler scheme. The vector potential is split into a source vector potential and a reduced vector potential which describes the eddy current and magnetization effects. The reduced formulation is compared with the underlying nonreduced formulation regarding the overall performance. In addition, the proper orthogonal decomposition method is applied to the iterative preconditioned conjugate gradient method to minimize the number of iterations per time step.


| INTRODUCTION
The finite element method (FEM) for magneto-quasistatic (MQS) field problems such as electric machines, transformers, and other power devices in the low frequency regime is commonly used for their design process. Typically, many MQS problems, also known as eddy current problems, involve materials with nonlinear magnetization characteristics. In order to solve transient eddy current problems implicit time integration schemes are widely used, since implicit time integration schemes are unconditionally stable. On the other hand, computationally expensive linearization schemes like the Newton-Raphson method are needed to carry out a single implicit time step. In Biro et al. 1 various FEM formulations for nonlinear transient eddy current problems using implicit time discretization, for example, the Crank-Nicolson scheme, are presented.
In order to avoid such linearization schemes, explicit time integration methods can be deployed. But unlike the implicit time integration, explicit time integration methods are stable by a maximum time step size, which must not be exceeded. In Yioultsis et al. 2 an explicit finite difference time domain approach is applied in the conductive regions, where the exterior nonconductive domain is computed using the boundary element method. Furthermore, in Ausserhofer et al. 3 a mixed time stepping procedure is proposed, where the discontinuous Galerkin method is used in conductive regions for a more efficient explicit time stepping. In the nonconductive region the problem is discretized in time implicitly in order to choose large time steps.
Based on a Schur complement reformulation of the eddy current problem in the entire computation domain, originally proposed by Clemens et al., 4 a differential algebraic equation (DAE) is transformed into a system of ordinary differential equations (ODE) in the conductive regions and into an algebraic equation in the nonconductive area. Explicit time integration of the ODE system leads to small maximum stable time step sizes. In Dutiné et al., [5][6][7] several techniques are proposed in order to optimize the computation of explicit time steps. Furthermore, in Kähne et al. 8 the explicit higher order Runge-Kutta-Chebyshev time integration method is used to increase the maximum stable time step size. A reduced vector potential formulation is for example used in Clemens and Weiland 9 where the vector potential is decomposed into a source vector potential and into a reduced vector potential, which describes the perturbation of the vacuum magnetic field by the presence of conductive and permeable materials. This formulation originates from Biro et al., 1 and is adapted to the finite integration technique. For implicit time stepping, it is shown, that the preconditioned conjugate gradient (PCG) method requires less iterations by solving for the reduced vector potential.
In this article, a similar reduced vector potential formulation is applied and discretized in time by the explicit Euler time integration scheme. The reduced vector potential formulation is compared with the nonreduced formulation presented in Dutiné et al. [5][6][7] by the numerical simulation of the TEAM 10 and the TEAM 21 a -2 eddy current problems.
The article is organized as follows. After this introduction, in Section 2 the corresponding mathematical formulation is presented. Section 3 shows the numerical results and a discussion about the limitations of the reduced and nonreduced vector formulation in practice.

| MATHEMATICAL FORMULATION
The magnetic flux density is described as B where κ is the electrical conductivity, ν is the reluctivity. The time dependent source current density is described as is the spatial current distribution and i S t ð Þ is the time dependent coil current. The reluctivity can be nonlinear in ferromagnetic materials and is then typically dependent on B ! . The corresponding variational formulation reads: find A ! such that The computational domain Ω decomposed in conductive Ω c and nonconductive Ω n subdomains with outer normal vector n ! .
holds true for all w where Ω is the spatial domain. The vector potential is discretized by a finite set of curl-conforming edge elements: where w Þand N denotes the total number of degrees of freedom (DOF) a i . The DOF assigned to the electrical conductive parts Ω c of the spatial domain Ω are stored in the coefficient vector a c whereas the coefficient vector a n holds the DOF assigned to nonconductive subdomain Ω n , see Figure 1. This reordering yields the DAE system with the regular conductivity mass matrix M c , the singular curl-curl reluctivity matrix K c a c ð Þ in conductive regions and the curl-curl matrix in the nonconductive regions K c . The matrix K cn couples the DOF of the conductive area with the neighboring DOF in the nonconductive area.
Applying the generalized Schur complement where K þ n is the Moore-Penrose pseudo-inverse matrix of the curl-curl matrix K n with nontrivial null space, Equation (4) decomposes in to an system of ODEs and an algebraic equation given by Rewriting Equation (6) yields and this system of ODEs can be integrated in time using the explicit Euler time integration method. The formulation above corresponds to the nonreduced vector potential A Ã formulation.

| The reduced vector potential formulation
The magnetic vector potential can be decomposed like where A ! s describes the source field excited by the source current density J ! and A ! s describes the magnetization and eddy current effects. 9 Using the approach in Equation (9), the nonreduced vector potential formulation Equation (8) can be rewritten as and can be resolved after Similar to Equation (8), Equation (11) can also be integrated in time with the explicit Euler time integration scheme. This formulation is dubbed as the reduced vector potential A r formulation.
The source vector potential a c,s can be determined a priori with the unified current distribution X ! S r ! in a single magneto-static computation or by using the law of Biot-Savart. Therefore, all matrix-vector multiplications involving a c,s in Equation (11) can be computed once in a single preprocessing step and are to be scaled only with the time dependent current i S t ð Þ. Since the solution vector a c varies in time the matrix K c a c ð Þ must be reassembled in each time step if the reluctivity is nonlinear.

| The TEAM 21 a -2 benchmark problem
The linear eddy current TEAM 21 a -2 benchmark problem proposed in Cheng et al., 10 is simulated in order to validate and compare the reduced and the nonreduced vector potential formulations. This problem contains two contrary exited coils with 300 turns each and the exciting currents Both coils are placed in front of a nonmagnetic and electrically conductive metal plate with the conductivity κ ¼ 1:3889 Á 10 6 S=m, the reluctivity is globally ν 0 , the reluctivity of vacuum, see Figure 2. The FEM simulations are executed on a mesh of about 540 000 tetrahedral elements. Using basis functions of second polynomial order results in 3.4 million DOF.
In order to perform a single explicit Euler time step for Equations (8) and (11), respectively, a repeated matrixvector product with the constant pseudoinverse K þ n must be evaluated. This forms a multiple right-hand side problem where j p represents K T cn a c or K T cn a c,r in Equations (6) and (11), respectively. This linear system of equations is solved using an iterative Jacobi-preconditioned CG method. 5,7 In order to solve for a single linear system according to Equation (13), the PCG method requires 970 iterations with a time consumption of 90 s for reduced vector potential formulations and 1400 iterations with a time consumption of 140 s for the nonreduced formulation. Without any further modifications, the reduced vector potential formulation provides over 30% of time saving over the nonreduced formulation. The explicit Euler time integration scheme is conditionally stable by a maximum time step size where λ max is the maximum eigenvalue of the matrix M À1 c K S a c ð Þ. For the spatial discretization of this problem Equation (14) leads to a very small maximum stable time step size Δt max . In order to simulate one sine period of 50 Hz according Equation (14) there are about 500 000 time steps to be solved. Apparently, both vector potential formulations are not practicable to perform a whole simulation in a reasonable amount of time.

| The TEAM 10 benchmark problem
The performance of the iterative PCG method can be enhanced by improved initial guess vectors. In this article, the proper orthogonal decomposition (POD) is used to compute initial values for the PCG method. 7,11,12 In addition, the nonlinear MQS TEAM 10 benchmark problem is simulated to compare both vector potential formulations in combination with the POD. 13 The problem geometry consists of symmetrically placed steel plates around a coil, see Figure 3. The steel plates have a nonlinear magnetization characteristic. The spatial mesh consists of 280 000 tetrahedral elements. The FEM simulation of 1st polynomial order has about 300 000 DOF. In Figure 4 is shown that both formulations provide almost identical results. Furthermore, using the POD, both the reduced and nonreduced formulation require about the identical CPU time to perform about 400 000 explicit Euler time steps in 150 min. The use of POD generated start vectors is thus shown to mostly eliminate the numerical advantage of the reduced vector potential formulation.
To investigate the performance of the POD in the nonreduced formulation more in detail, the time consumption and number of PCG iterations in each explicit Euler time step is recorded by the simulation of the TEAM 21 a -2 benchmark problem. Figure 5 shows very good agreement between numerical results and real measurements. In Figure 6 a histogram shows the probability of the PCG iterations of the entire simulation. Notable is, that the POD provides initial guesses which are accurate enough that the PCG method requires not a single iteration in most cases. Furthermore, Figure 7 shows the time consumption of all individual tasks in order to perform a single explicit Euler time step in the F I G U R E 2 Structure and mesh of the TEAM 21 a -2 benchmark problem. case of zero PCG iterations. Note that the PCG method still has the highest time consumption. Although there is no iteration performed the evaluation of the residual norm involves a computationally expensive matrix-vector multiplication.
All numerical simulations were executed on a workstation with an Intel XeonTM E5 processor using 16 parallel MPI processes. The FEM was implemented using the MFEM library. 14 The implementations for the PCG and POD methods were taken from the PETSc library. 15 F I G U R E 4 Numerical results of the reduced vector potential formulation (A r ) compared with the numerical results of the nonreduced vector potential formulation (A Ã ) for the TEAM 10 benchmark problem. Average magnetic flux computed at cross-sections S1-S3, see Figure 3.

| CONCLUSION
This article introduced a reduced vector potential formulation for transient eddy current problems where the reduced vector potential only describes magnetization and eddy current effects and the source vector potential only describes the air coil field. A generalized Schur complement transformed the spatially discretized DAE system into a system of ODEs which is integrated in time by the explicit Euler scheme to avoid the use of a linearization scheme as for example, the Newton-Raphson method.
The reduced vector potential formulation has shown a significant time saving over the nonreduced formulation by using the iterative PCG solver with the singular curl-curl matrix involved. Due to the massive oversampling related to the numerical stability limit of the explicit Euler time integration scheme the reduced formulation is still not practicable to perform a whole simulation in a reasonable amount of time.
It has been shown that improved initial guess vectors by the POD significantly reduced the number of required PCG iterations, often providing already sufficiently accurate approximative solution vectors resulting in zero PCG iterations. This can be considered for further investigations in order to optimize the total time consumption of this method. F I G U R E 6 Histogram for the PCG iterations in order to solve the linear system with the matrix K n in the TEAM 21 a -2 problem (second FE order) with improved initial guesses using the POD.
F I G U R E 7 Percentual time consumption of the individual tasks in order to perform a single explicit Euler time step in the A Ã formulation for the TEAM 21 a -2 problem (second FE order). Shown is the case if the PCG solver perform zero iterations in task 1) with POD enabled.