0 12pt 0 200pt 20pt Copyright © 2021 IEEE DOI: 10.1109/TMAG.2021.3057828

# A Two-Step Darwin Model Time Domain Formulation for Quasistatic
Electromagnetic Field Calculations
††thanks: Corresponding author: M. Clemens (email<EMAIL_ADDRESS>
This work was supported in parts by the Deutsche Forschungsgemeinschaft (DFG)
under grant no. CL143/11-2.

M. Clemens1, F. Kasolis1, M.-L. Henkel1, B. Kähne1, and M. Günther2  Institute
of Mathematical Modelling, Analysis and Computational Mathematics IMACM,
University of Wuppertal,  1 Chair of Electromagnetic Theory, Rainer-Gruenter-
Straße 21, 42119 Wuppertal, Germany  2 Chair of Numerical Analysis and Applied
Mathematics, Gaußstraße 20, 42119 Wuppertal, Germany

###### Abstract

In the absence of wave propagation, transient electromagnetic fields are
governed by a composite scalar/vector potential formulation for the
quasistatic Darwin field model. Darwin-type field models are capable of
capturing inductive, resistive, and capacitive effects. To avoid possibly non-
symmetric and ill-conditioned fully discrete monolithic formulations, here, a
Darwin field model is presented which results in a two-step algorithm, where
the discrete representations of the electric scalar potential and the magnetic
vector potential are computed consecutively. Numerical simulations show the
validity of the presented approach.

###### Index Terms:

Computational electromagnetics, electromagnetic fields, numerical simulation,
time domain analysis.

##  I Introduction

Electromagnetic field models that do not account for radiation effects are
dubbed quasistatic field models. For static fields, the Maxwell equations
decouple and enable the consideration of resistive, and capacitive or
inductive effects, with either electrostatic, or magnetostatic formulations,
separately. For capacitive-resistive effects, the electro-quasistatic (EQS)
field model is applicable, while resistive-inductive effects can be modelled
with the magneto-quasistatic (MQS) field approximation [1]. Quasistatic field
scenarios where inductive, resistive, and capacitive effects need to be
considered simultaneously, appear in high-frequency coils and coils of
inductive charging systems, where the capacitive effects between the coil
windings need to be taken into account and they are a common problem in
electromagnetic compatibility, e.g. in automotive engineering. For such
scenarios, it is quite common to use lumped $R$, $L$, $C$ parameter circuit-
type models such as Kirchhoff’s model or circuit models in combination with
field models used either for parameter extraction or in strong/weakly coupled
models, and circuit-formulation oriented partial-element equivalent circuit
(PEEC) methods. Rather recently, also field oriented models based on the
Darwin field formulation are considered. These quasistatic electromagnetic
field models are represented in terms of combined electric scalar and magnetic
vector potentials and feature a modified version of Ampère’s law by
eliminating the rotational parts of the displacement currents, i.e., by
neglecting the radiation effects in the model. Darwin formulations [2] are not
gauge-invariant, and thus, a number of different Darwin model formulations
have been considered, [3], [4], [5], [6], [7], [8], [9], [10], and [11]. The
paper is organized as follows. After this introduction, Darwin field models
with different established gauge conditions are highlighted. In the third
section, a Darwin model is presented which allows to use a two-step numerical
solution scheme. Section IV is comprised of numerical experiments with the
two-step time domain Darwin formulation, and is followed by conclusions.

##  II The Darwin Field Model

Darwin or Darwin-type field models for quasistatic electromagnetic field
distributions can be obtained by considering a decomposition of the electric
field intensity $\bm{E}$ into an irrotational part $\bm{E}_{\mathrm{irr}}$ and
a remaining part $\bm{E}_{\mathrm{rem}}$,

$\bm{E}=\bm{E}_{\mathrm{irr}}+\bm{E}_{\mathrm{rem}},$ (1)

where the irrotational part is represented as the gradient of an electric
scalar potential $\varphi$, that is,
$\bm{E}_{\mathrm{irr}}=-\operatorname{grad}\varphi$. The remainder part is
represented by the time derivative of a magnetic vector potential $\bm{A}$,
i.e., $\bm{E}_{\mathrm{rem}}=-\partial\bm{A}/\partial t$. Hence,

$\bm{E}=-\frac{\partial}{\partial
t}\bm{A}-\operatorname{grad}\varphi\text{,}\qquad\bm{B}=\operatorname{curl}\bm{A}$
(2)

holds for the electric field intensity and for the magnetic flux density,
respectively.

The assumption of a quasistatic electromagnetic field model enables the
elimination of the rotational parts of the displacement currents,
$\varepsilon\partial^{2}\bm{A}/\partial t^{2}\cong\bm{0}$, in Ampère’s law.
The result of this elimination is the so-called Darwin-Ampère equation

$\operatorname{curl}(\nu\operatorname{curl}\bm{A})+\kappa\frac{\partial}{\partial
t}\bm{A}+\kappa\operatorname{grad}\varphi+\varepsilon\operatorname{grad}\frac{\partial}{\partial
t}\varphi=\bm{J}_{\mathrm{S}},$ (3)

where $\nu$ is the reluctivity, $\kappa$ is the electric conductivity,
$\varepsilon$ is the permittivity, and $\bm{J}_{\mathrm{S}}$ is a source
current density.

The original formulation of the Darwin model [2] can be obtained by enforcing
a Coulomb gauge $\operatorname{div}\left(\bm{A}\right)=0$ corresponding to a
Helmholtz decomposition of the electric field intensity $\bm{E}$, i.e.,
assuming $\operatorname{curl}\bm{E}_{\mathrm{irr}}=0$ and
$\operatorname{div}\bm{E}_{\mathrm{rem}}=-\operatorname{div}\left(\partial\bm{A}/\partial
t\right)=0$. Intended to model charges in free space without conductive
materials, i.e., $\kappa=0$, $\varepsilon=\varepsilon_{0}$ and $\nu=\nu_{0}$,
as a consequence to the Helmholtz decomposition, the Gauß law does not
consider the rotational parts of the electric field and yields the
electrostatic Poisson equation as a gauge equation. Thus, the original Darwin
formulation is given by

$\displaystyle-\nu_{0}\Delta\bm{A}+\varepsilon\operatorname{grad}\frac{\partial}{\partial
t}\varphi$ $\displaystyle=\bm{J}_{\mathrm{S}},$ (4)
$\displaystyle\operatorname{div}\operatorname{grad}\varphi$
$\displaystyle=-\rho_{\mathrm{E}}/\varepsilon_{0},$ (5)

which requires to know the electric charge density $\rho_{\mathrm{E}}$ and its
motion with $\bm{J}_{\mathrm{S}}=\rho_{\mathrm{E}}\bm{v}$ along some velocity
vector $\bm{v}$.

To eliminate the free space assumption of the original Darwin model and to
include conductors, and permeable and dielectric materials, an application of
the divergence operator to the Darwin-Ampère equation (3) results in a
modified Darwin field formulation [5], [6], and yields the Darwin continuity
equation

$\operatorname{div}\left(\kappa\frac{\partial}{\partial
t}\bm{A}+\kappa\operatorname{grad}\varphi+\varepsilon\operatorname{grad}\frac{\partial}{\partial
t}\varphi\right)=\operatorname{div}\bm{J}_{\mathrm{S}}.$ (6)

This Darwin continuity equation lacks the radiation term
$\operatorname{div}(\varepsilon\partial^{2}\bm{A}/\partial t^{2})$, which is
present in the full Maxwell continuity equation. It was shown [5] that the
combined discrete formulation of the Darwin-Ampère equation (3) and the Darwin
continuity equation (6) results in non-symmetric systems. In addition, the
resulting system is singular and requires an additional gauge for the magnetic
vector potential in the non-conductive regions of the problem, such as
artificial conductivity [5] or an additional Coulomb-type gauge [8], [11]

$\operatorname{div}\left(\varepsilon\frac{\partial}{\partial
t}\bm{A}\right)=0,$ (7)

which can be enforced by adding this term with a scaling factor $1/\Delta t$
to the Darwin continuity equation (6), such that the gauge equation becomes a
temporally semi-discrete version of the full Maxwell continuity equation,
expressed in terms of the electrodynamic potentials $\bm{A}$ and $\varphi$.
The Coulomb-type gauge (7) can be additionally imposed as a third equation via
a Lagrange multiplier formulation [8]. Both Darwin model field formulations,
[8] and [11], are symmetric and do not require additional regularization.

##  III Two-Step Darwin Model Algorithms

The Darwin continuity equation (6) is extended with an additional gauge term
$\operatorname{div}\left(\kappa\partial\bm{A}/\partial t\right)=0$ [10] to
yield the electro-quasistatic equation

$\operatorname{div}\left(\kappa\operatorname{grad}\varphi+\varepsilon\operatorname{grad}\frac{\partial}{\partial
t}\varphi\right)=\operatorname{div}\bm{J}_{\mathrm{S}}.$ (8)

The expression $\operatorname{div}\left(\kappa\partial\bm{A}/\partial
t\right)$ omitted in (8) from the Darwin continuity equation (6), corresponds
to explicitly enforcing divergence-free eddy currents in conductive media,
i.e., neglecting eventually arising sources and sinks of current densities due
to the irrotational parts of the electric field. The combination of equation
(3) rewritten as

$\operatorname{curl}(\nu\operatorname{curl}\bm{A})+\kappa\frac{\partial}{\partial
t}\bm{A}=-\kappa\operatorname{grad}\varphi-\varepsilon\operatorname{grad}\frac{\partial}{\partial
t}\varphi+\bm{J}_{\mathrm{S}},$ (9)

with equation (8) results in a two-step formulation, where first the electro-
quasistatic total current density
$\bm{J}_{\mathrm{total}}=-\kappa\operatorname{grad}\varphi-\varepsilon\operatorname{grad}\left(\partial\varphi/\partial
t\right)+\bm{J}_{\mathrm{S}}$ is used as a solenoidal source term with
$\operatorname{div}\bm{J}_{\mathrm{Total}}=0$ to a magneto-quasistatic
formulation for the magnetic vector potential $\bm{A}$ represented by the
lefthand side of (9). This modified magneto-quasistatic formulation, however,
initially does not not address irrotational parts of $\bm{A}$ in the non-
conductive regions. While this does not affect the evaluation of $\bm{B}$ in
(2), the evaluation of the electric field according to (2) involves the
expression $\partial\bm{A}/\partial t$ also in the non-conductive regions,
which is commonly not covered in magneto-quasistatic field formulations.

To control the irrotational parts of $\bm{A}$, the magneto-quasistatic
formulation needs to be regularized. For this, the introduction of a small
artificial electrical conductivity $\hat{\kappa}$ in the non-conducting
regions has been suggested [5]. In case that $\kappa\gg 1/(\Delta
t)\varepsilon$ holds for a given time-step length $\Delta t$, a modified
electrical conductivity as e.g. $\hat{\kappa}=\kappa+1/(\Delta t)\varepsilon$
will regularize the formulation, where the resulting time-discrete
formulations will feature expressions of the type $1/(\Delta
t)\kappa+1/(\Delta t)^{2}\varepsilon$ as they occur in second-order time
discretization schemes, as e.g. Newmark-beta schemes used for full wave
Maxwell-Ampère equations [12]. Alternatively, a grad-div term augmentation for
spatially discretized magneto-quasistatic formulations is applicable [13],
[14].

The introduction of a small artificial electrical conductivity $\hat{\kappa}$
in the non-conducting regions is also an established technique to mitigate the
static limit instability of the electro-quasistatic formulation (8) that is
known to occur for $\frac{\partial}{\partial t}\varphi\rightarrow 0$.

By assuming that $\operatorname{div}\left(\kappa\partial\bm{A}/\partial
t\right)=0$ holds in equation (6), the calculation of the electric scalar
potential $\varphi$, using (8), is decoupled from that of the magnetic vector
potential $\bm{A}$. Thus, it is possible to independently first solve an
electro-quasistatic initial-boundary value problem that corresponds to (8),
and in a second step solve the modified magneto-quasistatic problem (9) with
the then available total current densities $\bm{J}_{\mathrm{Total}}$, as
depicted in Algorithm 1.

Algorithm 1 Two-Step Darwin Time Domain

1:Initialize $\varphi(t^{0})$ and $\bm{A}(t^{0})$;

2:for $n\leftarrow 0:n_{\mbox{\scriptsize End}}-1$ do

3: Solve problem (8) for $\varphi(t^{n+1})$;

4:end for

5:for $n\leftarrow 0:n_{\mbox{\scriptsize End}}-1$ do

6: Solve problem (9) for$\bm{A}(t^{n+1})$;

7: Evaluate $\bm{E}(t^{n+1})$ and $\bm{B}(t^{n+1})$ with (2);

8:end for

Alternatively, it is possible to consecutively execute a solution step for an
electro-quasistatic and a magneto-quasistatic field formulation for each
discrete timestep, using suitable time stepping schemes [15].

###  III-A Discrete Two-Step Darwin Time Domain Schemes

Reformulating (6) and (9) with a spatial volume discretization scheme, such as
the finite integration technique [16] or the finite element method with
Nédélec elements [17], results in a coupled system of time continuous matrix
equations

$\rm{\bf G}^{\\!\top}{\bf M}_{\kappa}\rm{\bf G}\bm{\upphi}+\rm{\bf
G}^{\\!\top}{\bf M}_{\varepsilon}\rm{\bf
G}\frac{\mathrm{d}}{\mathrm{d}t}\bm{\upphi}=\rm{\bf G}^{\\!\top}{\rm\bf
j}_{\mathrm{s}}-\rm{\bf G}^{\\!\top}{\bf
M}_{\kappa}\frac{\mathrm{d}}{\mathrm{d}t}{\rm\bf a},$ (10) $\mbox{\rm\bf
C}^{\\!\top}{\bf M}_{\nu}\mbox{\rm\bf C}{\rm\bf a}+{\bf
M}_{\kappa}\frac{\mathrm{d}}{\mathrm{d}t}{\rm\bf a}={\rm\bf
j}_{\mathrm{s}}-{\bf M}_{\kappa}\rm{\bf G}\bm{\upphi}-{\bf
M}_{\varepsilon}\rm{\bf G}\frac{\mathrm{d}}{\mathrm{d}t}\bm{\upphi},$ (11)

where ${\rm\bf a}$ is the degrees of freedom (dof) vector related to the
magnetic vector potential, $\bm{\upphi}$ is the dof vector of electric nodal
scalar potentials, ${\rm\bf j}_{\mathrm{s}}$ is a vector of transient source
currents, C is the discrete curl operator matrix, $\rm{\bf G}$ and $\rm{\bf
G}^{\\!\top}$ are discrete gradient and (negative) divergence operator
matrices. The matrices ${\bf M}_{\nu}$, ${\bf M}_{\kappa}$, ${\bf
M}_{\varepsilon}$ are discrete material matrices of possibly nonlinear
reluctivities, conductivities and permittivities, respectively, corresponding
to the specific discretization scheme in use. Employing e.g. an implicit Euler
backward differentiation time stepping scheme with time step $\Delta t$ to
(10) and (11) and ${\bf M}_{\sigma}={\bf M}_{\kappa}+(1/\Delta t){\bf
M}_{\varepsilon}$ yields a coupled system

$\displaystyle\left[\rm{\bf G}^{\\!\top}{\bf M}_{\sigma}\rm{\bf
G}\right]\bm{\upphi}^{n+1}$ $\displaystyle={\rm\bf f}_{1}({\rm\bf a}^{n+1}),$
(12) $\displaystyle\bigl{[}\mbox{\rm\bf C}^{\\!\top}{\bf M}_{\nu}\mbox{\rm\bf
C}+\frac{1}{\Delta t}{\bf M}_{\kappa}\bigr{]}{\rm\bf a}^{n+1}$
$\displaystyle={\rm\bf f}_{2}(\bm{\upphi}^{n+1}),$ (13)

where, using the notation $\Delta{\rm\bf a}^{n+1}={\rm\bf a}^{n+1}-{\rm\bf
a}^{n}$, the right-hand side vectors are

${\rm\bf f}_{1}={\rm{\bf G}}^{\\!\top}{\rm\bf j}_{s}^{n+1}+\frac{1}{\Delta
t}{\rm{\bf G}}^{\\!\top}{\bf M}_{\varepsilon}{\rm{\bf
G}}\bm{\upphi}^{n}-\frac{1}{\Delta t}{\rm{\bf G}}^{\\!\top}{\bf
M}_{\kappa}\Delta{\rm\bf a}^{n+1},$ (14) ${\rm\bf f}_{2}={\rm\bf
j}_{s}^{n+1}+\frac{1}{\Delta t}{\bf M}_{\kappa}{\rm\bf a}^{n}-{\bf
M}_{\sigma}{\rm{\bf G}}\bm{\upphi}^{n+1}+\frac{1}{\Delta t}{\bf
M}_{\varepsilon}{\rm{\bf G}}\bm{\upphi}^{n}.$ (15)

System (12), (13) are solved for each time step, starting from initial values
${\rm\bf a}^{0}={\rm\bf a}(t^{0})$ and $\bm{\upphi}^{0}=\bm{\upphi}(t^{0})$.
Adopting an iterative solution approach with iteration index $i=0,1,2,\ldots$
for each time step $t^{n+1}$ requires to provide an initial guess vector
${\rm\bf a}^{n+1}_{i=0}$ with ${\rm\bf f}_{1,i=0}={\rm\bf f}_{1}({\rm\bf
a}^{n+1}_{i=0})$. Inserting the solution vector of (12) rewritten as
$\bm{\upphi}^{n+1}=\left[\rm{\bf G}^{\\!\top}{\bf M}_{\sigma}\rm{\bf
G}\right]^{-1}{\rm\bf f}_{1,i=0}$ into the right-hand side vector equation
${\rm\bf f}_{2,i=0}={\rm\bf f}_{2}(\bm{\upphi}^{n+1}_{i=0})$ of (13) yields an
expression for the next iterative solution vector ${\rm\bf a}^{n+1}_{i=1}$.
Left application of the discrete divergence operator ${\rm{\bf G}}^{\\!\top}$
to this equation using the relation ${\rm{\bf G}}^{\\!\top}\mbox{\rm\bf
C}^{\\!\top}=\bm{0}$ yields the identity ${\rm{\bf G}}^{\\!\top}{\bf
M}_{\kappa}{\rm\bf a}^{n+1}_{i=1}={\rm{\bf G}}^{\\!\top}{\bf
M}_{\kappa}{\rm\bf a}^{n+1}_{i=0}$ and by induction

${\rm{\bf G}}^{\\!\top}{\bf M}_{\kappa}{\rm\bf a}^{n+1}_{i}={\rm{\bf
G}}^{\\!\top}{\bf M}_{\kappa}{\rm\bf a}^{n+1}_{0}\quad\forall
i\in\\{0,1,\ldots\\},$ (16)

i.e., in exact arithmetics the converged solution ${\rm\bf a}^{n+1}$ of the
iterative process will maintain the discrete divergence of its initial guess
solution ${\rm\bf a}^{n+1}_{i=0}$ in the right-hand side (14) of (12). An
initial guess ${\rm\bf a}^{n+1}_{i=0}={\rm\bf a}^{n}$ yields the result

${\rm{\bf G}}^{\\!\top}{\bf M}_{\kappa}{\rm\bf a}^{n}={\rm{\bf
G}}^{\\!\top}{\bf M}_{\kappa}{\rm\bf a}^{0}\quad\forall n\in\\{0,1,\ldots\\}.$
(17)

Thus, with the choice of the initialization vector ${\rm\bf a}^{0}$ of the
time integration process at $t^{0}$ acting as an initial gauge, the difference
expression ${\rm{\bf G}}^{\\!\top}{\bf M}_{\kappa}\bigl{[}{\rm\bf
a}^{n+1}-{\rm\bf a}^{n}\bigr{]}=\bm{0}$ in (14) vanishes for all time steps,
and thus, system (12), (13) gets decoupled.

##  IV Numerical Experiments

To verify the performance of the proposed two-step Darwin time domain
algorithm, two three-dimensional copper coil problems are considered
(electrical conductivity $\kappa_{\mathrm{copper}},=5.96\cdot
10^{7}~{}\mathrm{S/m}$. Both problems are illustrated in Fig. 1. For each
case-study, the computational domain is
$\Omega=(\Omega_{0}\cup\overline{\Omega}_{\kappa})\setminus(\Gamma_{\mathrm{E}}\cup\Gamma_{\mathrm{G}})\subset\mathbb{R}^{3}$
and is free from charge and current sources. The bounding surfaces are
perfectly conducting, with $\Gamma_{\mathrm{G}}$ being grounded and
$\Gamma_{\mathrm{E}}$ supplying the transient excitation

$\varphi(t)=\varphi_{\mathrm{max}}\cdot f\cdot\min(t,1/f)\cdot\sin(2\pi ft),$
(18)

where $\varphi_{\mathrm{max}}=12~{}\mathrm{V}$ is the maximum voltage and
$f=10~{}\mathrm{MHz}$ is the excitation frequency with a wavelength
$\lambda=30~{}\mathrm{m}$ in void. The longest side of the domain $\Omega$
associated with the helical coil is $l=6.3~{}\mathrm{cm}$, the one associated
with the planar coil is $l=1.35~{}\mathrm{cm}$. Since $\ell\ll\lambda$ in both
cases, the radiation-free assumption of the Darwin field model is justified.

Figure 1: In both coil case-studies, $\Omega_{0}$ is void, with $\hat{\kappa}$
being a small artificial electrical conductivity ($\hat{\kappa}=10^{-2}$ S/m),
while $\Omega_{\kappa}$ is occupied by copper.

The problems that constitute the two-step algorithm are discretized in space
with the FEM, using first-order Lagrange elements for the scalar electric
problem and zeroth-order Nédélec elements for the vectorial magnetic problem;
see Table I for the number of degrees of freedom in each finite element space.
Regarding time-discretization, both problems have been integrated with the
trapezoidal method, which is implicit, second-order accurate, and A-stable,
with different time steps $\Delta t\in\\{2.5,1.25,0.625\\}~{}\mathrm{ns}$ for
a total of $t_{\mathrm{End}}=(n_{\mathrm{End}}-1)\Delta t=1200~{}\mathrm{ns}$.
With the excitation functions in (18), the two-step algorithm is expected to
yield an approximation of a frequency-domain full Maxwell solution, and hence,
the latter is used for obtaining reference solutions on the same meshes. For
all linear systems a direct solver is used.

TABLE I: Number of Degrees of Freedom (dof) | Lagrange Elements | Nédélec Elements
---|---|---
Helical Coil | $16\,882$ | $118\,609$
Planar Coil | $41\,357$ | $292\,905$

In Fig. 2, the magnetic flux density and the electric field intensity are
depicted for each coil. There, the two-step algorithm for the Darwin field
model, successfully captures, not only the induction, but also the capacitance
between the coil windings. Table II depicts the relative differences of the
electric field approximations provided by the Darwin field model for the
planar coil model at frequencies ranging from $10~{}\mathrm{kHz}$ up through
$1~{}\mathrm{GHz}$. The results show that the remainder part
$-\partial\bm{A}/\partial t$ of the electric field needs to be evaluated in
(2) which necessitates the regularization of (9).

Figure 2: The magnitude of the magnetic flux density and the electric field
intensity for the two coil problems at $t=825~{}\mathrm{ns}$.

In Fig. 3, the difference between the Maxwell and Darwin field models is
quantified with the norm

$\|\operatorname{Re}(\bm{F}_{\mathrm{M}})-\bm{F}_{\mathrm{D}}\|_{L^{2}(\Omega)}/\|\bm{F}_{\mathrm{M}}\|_{L^{2}(\Omega)},$
(19)

where $\bm{F}\in\\{\bm{B},\bm{E}\\}$ is a physical field quantity, computed as
in (2), and the subscripts $\mathrm{M}$ and $\mathrm{D}$ stand for the Maxwell
and Darwin field models, respectively. In Fig. 3, the first row of results is
associated with the helical coil, while the second row with the planar coil.
In the same figure, the effect of the time discretization scheme is also
apparent, with a tendency towards improved accuracy for smaller time steps,
since convergence to the time-harmonic solution is expected.

Figure 3: The relative difference between the field quantities $\bm{B}$ and $\bm{E}$ for the helical coil (first row) and for the planar coil (second row), computed using a full Maxwell frequency-domain solver and the two-step Darwin algorithm. TABLE II: Relative $\bm{E}$-field Differences at $t=3.125/f$ for the Planar Coil as Functions of the Frequency $f$. $f$ [Hz] | $\|\operatorname{Re}(\bm{E}_{\mathrm{M}})-\bm{E}_{\mathrm{D}}\|/\|\bm{E}_{\mathrm{M}}\|$ | $\|\operatorname{Re}(\bm{E}_{\mathrm{M}})-\bm{E}_{\mathrm{D,irr}}\|/\|\bm{E}_{\mathrm{M}}\|$
---|---|---
$10^{4}$ | $8.22\cdot 10^{-6}$ | $2.74\cdot 10^{-3}$
$10^{5}$ | $9.31\cdot 10^{-6}$ | $2.92\cdot 10^{-2}$
$10^{6}$ | $4.60\cdot 10^{-5}$ | $1.59\cdot 10^{-1}$
$10^{7}$ | $7.33\cdot 10^{-4}$ | 1.47
$10^{8}$ | $2.00\cdot 10^{-3}$ | 2.28
$10^{9}$ | $1.41\cdot 10^{-2}$ | 2.22

##  V Conclusions

A two-step algorithm for the transient $(\bm{A},\varphi)$ formulation of the
quasistatic Darwin field model is introduced and, for the first time,
numerically validated against the full system of Maxwell’s equations. The
presented two-step Darwin time domain quasistatic field formulation accounts
for capacitive, inductive, and resistive effects. The advantages of this
scheme result from consecutively combining an electro-quasistatic and a
modified magneto-quasistatic field model, and thus, it benefits from existing
efficient time domain solution techniques that provide flexibility in terms of
material non-linearities and excitation profiles.

## References

  * [1] H. A. Haus and J. R. Melcher, Electromagnetic Fields and Energy.  Englewood Cliffs, New Jersey: Prentice Hall, 1989.
  * [2] C. G. Darwin, “The dynamical motion of particles,” Phil. Mag., vol. 93, pp. 537 – 551, 1920.
  * [3] P. A. Raviart and E. Sonnendrücker, “Approximate models for the Maxwell equations,” J. Comp. Appl. Math., vol. 63, pp. 69 – 81, 1995.
  * [4] J. Larsson, “Electromagnetics from a quasistatic perspective,” Am. J. Phys., vol. 75, no. 3, pp. 230–239, 1995.
  * [5] S. Koch and T. Weiland, “Different types of quasistationary formulations for time domain simulations,” Radio Science, vol. 46, 2011.
  * [6] S. Koch, H. Schneider, and T. Weiland, “A low-frequency approximation to the Maxwell equations simultaneously considering inductive and capacitive phenomena,” IEEE Trans. Magn., vol. 48, no. 2, pp. 511–514, 2012.
  * [7] I. C. Garcia, S. Schöps, H. D. Gersem, and S. Baumanns, “Chapter 1 systems of differential algebraic equations in computational electromagnetics,” in Applications of Differential-Algebraic Equations: Examples and Benchmarks (S. Campell, A. Ilchmann, V. Mehrmann, and T. Reis, eds.), pp. 123–169, Springer Verlag, 2018.
  * [8] Y. Zhao and Z. Tang, “A novel gauged potential formulation for 3d electromagnetic field analysis including both inductive and capacitive effects,” IEEE Trans Magn., vol. 55, pp. 1–5, June 2019.
  * [9] Z. Badics, S. Bilicz, J. Pávo, and S. Gyimóthy, “Finite element $a-v$ formulations for quasistatic Darwin models,” in Proc. IEEE CEFC 2018 Conference, Hangzhou, China, 2018.
  * [10] M. Clemens, B. Kähne, and S. Schöps, “A Darwin time domain scheme for the simulation of transient quasistatic electromagnetic fields including resistive, capacitive and inductive effects,” in 2019 Kleinheubach Conference, Miltenberg, Germany, pp. 1–4, 2019.
  * [11] H. Kaimori, T. Mufune, and A. Kameari, “Novel application of Coulomb gauge condition in electromagnetic fem computations for Darwin approximation,” in IEEE CEFC 2020 Conference, Pisa, Italy, p. 497, 2020.
  * [12] D. C. Dibben and R. Metaxas, “Frequency domain vs. time-domain finite element methods for calculation of fields in multimode cavities,” IEEE Trans. Magn., vol. 33, no. 2, pp. 1468–1471, 1997.
  * [13] A. Bossavit, “’Stiff’ problems in eddy-current theory and the regularization of Maxwell’s equations,” IEEE Trans. Magn., vol. 37, pp. 3542–3545, Sept. 2001.
  * [14] M. Clemens, S. Schöps, H. De Gersem, and A. Bartel, “Decomposition and regularization of nonlinear anisotropic curl-curl daes,” COMPEL, vol. Vol. 30, no. No. 6, pp. 1701–1714, 2011.
  * [15] M. Clemens, “Large systems of equations in a discrete electromagnetism: Formulations and numerical algorithms,” IEE Proc. - Scie. Meas. Technol., vol. 152, no. 2, pp. 50–72, 2005.
  * [16] T. Weiland, “Time domain electromagnetic field computation with finite difference methods,” Int. J. Num. Mod.: ENDF, vol. 9, pp. 259–319, 1996\.
  * [17] J. C. Nédélec, “Mixed finite elements in ${R}^{3}$,” Numer. Math., vol. 35, pp. 315–341, 1980.