A Darwin Time Domain Scheme for the Simulation of Transient Quasistatic Electromagnetic Fields Including Resistive, Capacitive and Inductive Effects
cc (cid:13) A Darwin Time Domain Scheme for the Simulationof Transient Quasistatic Electromagnetic FieldsIncluding Resistive, Capacitive and Inductive Effects st Markus Clemens
University of WuppertalChair of Electromagnetic Theory
Wuppertal, [email protected] nd Bernhard K¨ahne
University of WuppertalChair of Electromagnetic Theory
Wuppertal, [email protected] rd Sebastian Sch¨ops
Technische Universit¨at DarmstadtGraduate School of Computational Engineering
Darmstadt, [email protected]
Abstract —The Darwin field model addresses an approximationto Maxwell’s equations where radiation effects are neglected.It allows to describe general quasistatic electromagnetic fieldphenomena including inductive, resistive and capacitive effects.A Darwin formulation based on the Darwin-Amp`ere equationand the implicitly included Darwin-continuity equation yieldsa non-symmetric and ill-conditioned algebraic systems of equa-tions received from applying a geometric spatial discretizationscheme and the implicit backward differentiation time integrationmethod. A two-step solution scheme is presented where theunderlying block-Gauss-Seidel method is shown to change theinitially chosen gauge condition and the resulting scheme onlyrequires to solve a weakly coupled electro-quasistatic and amagneto-quasistatic discrete field formulation consecutively ineach time step. Results of numerical test problems validate thechosen approach.
Index Terms —Electromagnetic fields, linear algebra, numericalsimulation, time domain analysis.
I. I
NTRODUCTION
Quasistatic field models derived from Maxwell’s equationsare considered valid, if the shortest wavelength of a problemwell exceeds the diameter of the considered problem [1], [2].In these cases, radiation effects can be neglected. For furthertaxonomy, the electric energy density w e and magnetic energydensity w m are considered: in case of w e (cid:29) w m everywhere inthe problem domain, the electro-quasistatic model is applicableand a variation of the magnetic electric can be neglected, i.e., ∂∂t B ≈ . Thus, the electric field is irrotational (curl E = ) and governed by resistive and capacitive effects. In case of w e (cid:28) w m , the magneto-quasistatic field approximation onlytakes into account resistive and inductive field effects wheredisplacement currents are neglected within Amp`ere’s law, i.e., ∂∂t D ≈ . Quasistatic field scenarios, where w e ≈ w m holds, i.e.,where capacitive, resistive and inductive field effects are tobe considered in the same problem, are often described usingthe full set Maxwell’s equations. As a result, in these modelsthe otherwise negligible radiation effects are still considered asan unnecessary part of the model. Especially in time domainformulations this results in a high stiffness of the resulting discrete field formulations. Alternatively, quasistatic modelsof such scenarios often involve the use of lumped parameterformulations based on Kirchhoff’s equations.The Darwin field model is an approximation to Maxwell’sequations related to general quasistatic field scenarios includ-ing capacitive, resistive and inductive field effects, i.e., whereonly radiation effects can be neglected [3], [4], [5], [6], [7],[8], [9].Following the notation in [8], in the quasistatic Darwin fieldmodel, the electric field E is subject to a decomposition E = E irr. + E rem. , (1)and split up into an irrotational part E irr. with curl( E irr. ) = , for which a scalar electric potential representation E irr. = − grad( ϕ ) exists, and a remainder part E rem. In (1), div( E rem. ) = only holds in the uniquely special case of aHelmholtz decomposition that is typically only admissible forhomogeneous material distributions. Alternatively, div( E rem. )may be nonzero, which needs to be taken into account withinan additional gauge. Within the Darwin field model, the rota-tional parts of the displacement current densities are neglectedwith ∂∂t ( ε E rot ) = . This essentially eliminates the hyperboliccharacter of the full Maxwell’s equations which is responsiblefor modeling wave propagation phenomena and translatesthe Darwin model into a system with only first order timederivatives. Based on these equations several reformulations interms of the magnetic vector potential A and scalar potential ϕ can be formulated.Following this introduction, in section II a quasistatic fieldmodel for time domain problems is formulated featuring theDarwin-Amp`ere’s equation and its corresponding Darwin-continuity equation in terms of electrodynamic potentials.Section III describes the space and time discretzation ofthis Darwin formulation resulting in ill-conditioned and non-symmetric monolithic algebraic systems of equations. SectionIV introduces a two-step solution technique which requiresonly symmetric algebraic systems to be solved. In Section V, a r X i v : . [ c s . C E ] M a r umerical test results are shown, including a discussion of theresults, followed by a conclusion.II. ( A , ϕ ) F ORMULATION FOR THE D ARWIN M ODEL
The splitting of the electric field in (1) is expressed in termsof electrodynamic potentials, i.e., the magnetic vector potential A and the scalar electric potential ϕ , with E rem. = − ∂∂ t A and E irr. = − grad ϕ. Initially, the electric field then reads E = − ∂∂ t A − grad ϕ. (2)The Darwin formulation ignores the rotational parts ofthe displacement current densities related to the radiationof electromagnetic waves, i.e., in the electric displacementcurrents , i.e., ∂ ∂t A ≈ . With this ansatz, the Amp`ere’s equation reduces to theDarwin-Amp`ere’s equation curl( ν curl A )+ κ ∂∂ t A + κ grad ϕ + ε grad ∂∂ t ϕ = J s , (3)where ν is the reluctivity, κ the specific electric conductivity, ε the permittivity and J s denotes a transient source currentdensity.For the solution of (3) an additional equation is required todescribe the relation of the magnetic vector potential A andthe electric scalar potential ϕ. Due to the gauge invariance ofthe Darwin field model, various Darwin formulations can bederived [4], [8], [9]. If non-homogeneous material distributionsneed to be considered in the field model, the coupling of A and ϕ in (3) by use of the implicitly contained continuityequation appears as an obvious choice. Left application of thedivergence operator to the Darwin-Amp`ere equation (3) yieldsa Darwin-model continuity equation [7]div (cid:18) κ ∂∂ t A + κ grad ϕ + ε grad ∂∂ t ϕ (cid:19) = div J s . (4)It should be noted, that (4) is not identical to the originalcontinuity equation div J + ∂∂ t ρ = related to the fullMaxwell model, where ρ is the electric space charge. Includingthe Gauß’ law, the full Maxwell model continuity equationcontains the expression ε ∂ ∂t ( A ) related to the rotational partsof the displacement currents. These are specifically omittedwithin the Darwin model.III. A D ISCRETE D ARWIN M ODEL F ORMULATION
For the general case of non-homogeneous material distri-butions, a spatial discretization to the equations (3) and (4)is considered following [7]. The application of a mimeticdiscretization scheme as e.g. the finite integration technique(FIT) [10], [11], Whitney finite element method (WFEM) [12]or the cell method [13] yields the systems of matrix equations C (cid:62) M ν C a + M κ dd t a + M κ G ϕ + M ε G ddt ϕ = j s , (5) G (cid:62) M κ ddt a + G (cid:62) M κ G ϕ + G (cid:62) M ε G ddt ϕ = G (cid:62) j s , (6)where a is the degrees of freedom (dof) vector related to themagnetic vector potential, ϕ is the dof vector of electric nodal scalar potentials, C is the discrete curl operator matrix, G and G (cid:62) are discrete gradient and (negative) divergence operatormatrices. The matrices M ν , M κ , M ε are the (possibly nonlin-ear) discrete material matrices of reluctivities, conductivitiesand permittivities, respectively, and the construction of thesediscrete Hodge operators depends on the specific discretizationscheme.The discrete Darwin equations (5) and (6) can be rewrittenas a first order differential-algebraic system of equations (cid:20) M κ M ε GG (cid:62) M κ G (cid:62) M ε G (cid:21) dd t (cid:20) a ϕ (cid:21) + (cid:20) C (cid:62) M ν C M κ G G (cid:62) M κ G (cid:21) (cid:20) a ϕ (cid:21) = (cid:20) j s G (cid:62) j s (cid:21) . (7)While the non-symmetry of the block matrices in (7) canbe partially eliminated, the presence of metallic objects in theproblem domain will result in large differences in the order ofmagnitude of the entries.A time discrete Darwin model based on the monolithic timedomain formulation (7) results from the application of e.g. afirst order convergent Euler backward differentiation (BDF1)formula [14] with a unconditionally stable time step ∆ t := t n +1 − t n . With M σ := M κ + t M ε the implicit time steppingscheme requires to solve the algebraic system of equations (cid:20) C (cid:62) M ν C + t M κ M σ G t G (cid:62) M κ G (cid:62) M σ G (cid:21) (cid:20) a ϕ (cid:21) n +1 = 1∆ t (cid:20) M κ M ε GG (cid:62) M κ G (cid:62) M ε G (cid:21) (cid:20) a ϕ (cid:21) n + (cid:20) j s G (cid:62) j s (cid:21) n +1 (8)for each time step and Newton iteration in case on nonlinearmaterial laws. The real-valued system matrix is non-symmetricand can not be symmetrized; it is singular, if M κ is singular.In [14] the system (8) is additionally regularized assuming anon-physical conductivity in the non-conductive regions of thesimulation problem. The system matrices of the monolithictime discrete algebraic systems of equations (8) can be ex-tremely ill-conditioned due to their off-diagonal matrix blockswith entries varying by different orders of magnitude. In [14] adirect solver based on sparse LU-decomposition was used fora discrete problem with a total number of only 12,101 dofs.IV. A T WO -S TEP D ARWIN T IME D OMAIN S CHEME
Rewriting the mutually coupled equations (5) and (6) into C (cid:62) M ν C a + M κ dd t a = j s − M κ G ϕ − M ε G ddt ϕ , (9) G (cid:62) M κ G ϕ + G (cid:62) M ε G ddt ϕ = G (cid:62) j s − G (cid:62) M κ ddt a (10)shows both the discrete magneto-quasistatic ( A (cid:63) ) formula-tion (9) (see also [15], [16], [17]) and the discrete electro-quasistatic scalar electric potential formulation (10) (see [18])coupled to each other with their specific right hand sidevectors, respectively.This motivates adopting a strong coupled iteration approachfor the solution of the time and space discretized reformu-lations of (9) and (10). This approach is identical to anterative block-Gauss-Seidel solution of the monolithic system(8) whose convergence can be shown by inspecting the spectralproperties of the matrices.Let’s denote by a n +1 m and ϕ n +1 m the iterative solution afterthe m -th Gauss-Seidel iteration at time step n +1 . The iterationscheme requires an initial guess, e.g. by extrapolation a n +10 := a n , (11)where a n denotes the final solution at time step n , i.e., after m = 1 , , . . . iterations.Due to discrete conservation properties of the discrete fieldformulations [11], the discrete divergence relation G (cid:62) M κ a n+1m+1 = G (cid:62) M κ a n+1m (12)holds. When using (11), this discrete conservation property(12) implicitly eliminates the time discrete backward deriva-tive expression t G (cid:62) M κ (cid:2) a n+1m+1 − a n (cid:3) within the iterationscheme. This translates to div (cid:0) κ ∂∂ t A (cid:1) = 0 in the continuouscase, i.e., the rotational parts of the eddy current densitieswill be solenoidal. This is an physically acceptable modelassumption in case of highly conductive materials withinwhich the effects of the displacement currents are negligible.Due to this, the Darwin continuity equation (4) reduces to thestandard electro-quasistatic field formulation [18].This property weakens the coupling of (9) and (10) and theiterative scheme reduces to the two-step approach shown inAlgorithm 1, where in each time step only two symmetric andpositive (semi-)definite systems of algebraic equations needto be solved consecutively. For this task efficient solutionschemes are available (see e.g. [17], [18], [19]). Algorithm 1
Two-Step Darwin Time Domain Algorithm Initialize ϕ := ϕ ( t ); a := a ( t ); for n ← n End : do Solve: (cid:2) G (cid:62) M σ G (cid:3) ϕ n +1 = G (cid:62) j n +1 s + t G (cid:62) M ε G ϕ n ; Solve: (cid:2) C (cid:62) M ν C + t M κ (cid:3) a n +1 = j n +1 s + t M κ a n − M σ G ϕ n +1 + t M ε G ϕ n ; end for V. N
UMERICAL R ESULTS
Two test structures are excited at f = 100 MHz, withdimensions small against the wave length, i.e., the quasistaticassumption holds. Using ramped sinusoidal excitations, thetwo-step Darwin time domain scheme is implemented usingthe MFEM library [21]. For the solution of the electro-quasistatic system and the weakly coupled magneto-quasistaticsystems at each time step, efficient algebraic multigrid (AMG)schemes provided in the PETSc [19] linear algebra solverlibrary are used.
A. High-frequency coil
A high-frequency coil structure of 63 mm length (see Fig. 1(left) is considered using a mesh consisting of 237,835 tetra-hedra. For the time-domain simulation a ramped sinusoidalexcitation profile is used for 10 periods. The dof vector ϕ hasdimension 39,038 and the vector a has dimension , .Fig 2 shows results achieved for the magnetic and the electricfield after t = 91 /f and t = 9 , /f . Fig. 1. Coil geometry (left); RLC structure geometry (right).Fig. 2. Magnitudes of the electric field E ( t = 9 T ) (top left) and magneticfield B ( t = 9 T ) (bottom left) and the electric field E ( t = 9 , T ) (topright) and magnetic field B ( t = 9 . T ) (bottom right) simulated with thetwo-step Darwin time domain algorithm. B. RLC Model
A RLC-structure of 63 mm length presented in [20], [9],i.e., a wire (electrical conductivity κ = 10 S/m connectinga coil and a capacitor with a dielectric inset ( ε r = 2 ) (seeFig. 1 (right)) is considered. The FEM mesh consisting of543,783 tetrahedra. For the time-domain simulation a rampedsinusoidal excitation profile is used for 10 periods. The dofvector ϕ has dimension , and the vector a has dimension , . Fig 3 shows the simulation results achieved with thewo-step Darwin time domain scheme for the magnetic andthe electric field. Fig. 3. Magnitudes of the electric field E ( t = 9 T ) (top left) and magneticfield B ( t = 9 T ) (bottom left) and the electric field E ( t = 9 , T ) (topright) and magnetic field B ( t = 9 . T ) (bottom right) simulated with thetwo-step Darwin time domain algorithm. C. Discussion
The simulation results achieved with the two-step timedomain method, show that resistive, inductive and capacitiveeffects are included: The capacitive coupling of the high-frequency coil windings in the HF coil is included with theirrotational parts of the electric field and shown in Fig. 2. Theresults in Fig. 2 and Fig. 3 show the exchange of electricand magnetic field energy in the sinusoidally excited teststructures. VI. C
ONCLUSION
The Darwin field model was analyzed to describe generalquasistatic electric and magnetic field distributions by onlyneglecting the rotational contributions of the displacementcurrents. Starting from an ( A , ϕ ) formulation of the Darwin-Amp`ere law and the Darwin-continuity equation, the resultingdiscrete Darwin model represents a differential-algebraic setof equations. In order to avoid the solution of ill-conditionedand non-symmetric monolithic algebraic systems of equationsrequired within implicit time discretization schemes, a two-step solution schemes was presented based on the consecutivesolution of weakly coupled discrete electro- and magneto-quasistatic field formulations in each time step. Numericalresults of quasistatic electromagnetic structures, where capaci-tive, inductive and resistive field effects need to be considered,showed the validity of this two-step approach and its abilityto solve realistic 3D problem resolutions by enabling the useof efficient solution schemes. A CKNOWLEDGEMENT
This work is supported in parts by the Deutsche Forschungs-gemeinschaft (DFG) under grant no. CLE143/10-2.R
EFERENCES[1] T. Steinmetz, S. Kurz, and M. Clemens, “Domains of validity ofquasistatic and quasistationary field approximations,”
COMPEL , vol. 30,no. 5, pp. 1237–1247, 2011.[2] V. Mazauric, L. Rondot, and P. Wendling, “Enhancing quasi-staticmodeling: A claim for electric field computation,”
IEEE Trans. Magn. ,vol. 49, pp. 1629–1632, May 2013.[3] P. A. Raviart and E. Sonnendr¨ucker, “A hierarchy of approximate modelsfor the Maxwell equations,”
Num. Math. , vol. 73, pp. 329 – 372,February 1996.[4] J. Larsson, “Electromagnetics from a quasistatic perspective,”
Am. J.Phys. , vol. 75, no. 3, pp. 230–239, 1995.[5] C. Liao and L. Ying, “An analysis of the Darwin model of approxi-mations to Maxwell’s equations in 3-d unbounded domains,”
Commun.Math. Sci. , vol. 6, no. 3, pp. 695–710, 2008.[6] N. Fang, C. Liao, and L. Ying, “Darwin approximation to Maxwell’sequations,” No. I, 2009. ICCS 2009.[7] S. Koch, H. Schneider, and T. Weiland, “A low-frequency approximationto the Maxwell equations simultaneously considering inductive andcapacitive phenomena,”
IEEE Trans. Magn. , vol. 48, pp. 511–514,February 2012.[8] I. Cortes-Garcia, S. Sch¨ops, H. D. Gersem, and S. Baumanns, “Chapter 1systems of differential algebraic equations in computational electromag-netics,” in
Applications of Differential-Algebraic Equations: Examplesand Benchmarks (S. Campell, A. Ilchmann, V. Mehrmann, and T. Reis,eds.), pp. 123–169, Springer Verlag, 2018.[9] Z. Badics, S. Bilicz, J. P´avo, and S. Gyim´othy, “Finite element a − v formulations for quasistatic darwin models,” in Proc. CEFC2018 Con-ference, Hangzhou, China , 2018.[10] T. Weiland, “A discretization method for the solution of Maxwell’sequations for six-component fields,”
Electronics and CommunicationsAE ¨U , vol. 31, no. 3, pp. 116–120, 1977.[11] M. Clemens and T. Weiland, “Discrete electromagnetism with theFinite Integration Technique,” in
Geometric Methods for ComputationalElectromagnetics (F. L. Teixeira, ed.), no. 32 in PIER, pp. 65–87,Cambridge, Massachusetts, USA: EMW Publishing, 2001.[12] J. C. N´ed´elec, “Mixed finite elements in R ,” Numer. Math. , vol. 35,pp. 315–341, 1980.[13] E. Tonti, “Finite formulation of the electromagnetic field,” in
GeometricMethods for Computational Electromagnetics (F. L. Teixeira, ed.), no. 32in Progress in Electromagnetic Research (PIER), pp. 1–44, Cambridge,Massachusetts, USA: EMW Publishing, 2001.[14] S. Koch and T. Weiland, “Different types of quasistationary formulationsfor time domain simulations,”
Radio Science , vol. 46, no. RS0E12, 2011.[15] C. R. I. Emson and J. Simkin, “An optimal method for 3-D eddycurrents,”
IEEE Trans. Magn. , vol. Mag-19, no. 6, pp. 2450–2452, 1983.[16] A. Kameari, “Calculation of transient 3d eddy current using edgeelements,”
IEEE Trans. Magn. , vol. 26, no. 5, pp. 466–469, 1990.[17] M. Clemens and T. Weiland, “Transient eddy current calculation with theFI-method,”
IEEE Trans. Magn. , vol. 35, no. 3, pp. 1163–1166, 1999.[18] M. Clemens, M. Wilke, G. Benderskaya, H. De Gersem, W. Koch, andT. Weiland, “Transient electro-quasi-static adaptive simulation schemes,”
IEEE Trans. Magn. , vol. 240, no. 2, pp. 1294–1297, 2004.[19] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, “Efficienctmanagement of parallelism in object oriented numerical software li-braries,” in
Modern Software Tools in Scientific Computing (E. Arge,A. M. Bruaset, and H. P. Langtangen, eds.), p. 163–202, BirkhauserPress, 1997.[20] M. Jochum, O. Farle, and R. Dyczij-Edlinger, “A new low-frequencystable potential formulation for the finite-element simulation of emfields,”