# Numerical Simulation-Based Investigation of the Limits of Different Quasistatic Models

Houssein Taha, Zuqi Tang, Thomas Henneron, Yvonnick Le Menach,

Florentin Salomez, Jean-Pierre Ducreux

► To cite this version:

Houssein Taha, Zuqi Tang, Thomas Henneron, Yvonnick Le Menach, Florentin Salomez, et al.. Numerical Simulation-Based Investigation of the Limits of Different Quasistatic Models. Applied Sciences, 2021, Applied Sciences, 11 (23), pp.11218. 10.3390/app112311218 . hal-03959278

HAL Id: hal-03959278

https://hal.univ-lille.fr/hal-03959278

Submitted on 27 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Numerical Simulation-Based Investigation of the Limits of Different Quasistatic Models

Houssein Taha<fnref n="1" />, Zuqi Tang<fnref n="1" />,* D, Thomas Henneron<fnref n="1" />, Yvonnick Le Menach<fnref n="1" />, Florentin Salomez<fnref n="1" /> and Jean-Pierre Ducreux<fnref n="2" />

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Taha, H.; Tang, Z.; Henneron, T.; Le Menach, Y.; Salomez, F.; Ducreux, J.-P. Numerical Simulation-Based Investigation of the Limits of Different Quasistatic Models. *Appl. Sci.* **2021**, *11*, 11218.
https://doi.org/10.3390/app112311218

Academic Editors: David Pardo and
Luis E. García-Castillo

Received: 30 September 2021

Accepted: 21 November 2021

Published: 25 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** The modeling of the capacitive phenomena, including the inductive effects becomes critical, especially in the case of a power converter with high switching frequencies, supplying an electrical device. At a low frequency, the electro-quasistatic (EQS) model is widely used to study the coupled resistive-capacitive effects, while the magneto-quasistatic (MQS) model is used to describe the coupled resistive-inductive effects. When the frequency increases, the Darwin model is preferred, which is able to capture the coupled resistive-capacitive-inductive effects by neglecting the radiation effects. In this work, we are interested in specifying the limits of these models, by investigating the influence of the frequency on the electromagnetic field distributions and the impedance of electromagnetic devices. Two different examples are carried out. For the first one, to validate the Darwin model, the measurement results are provided for comparison with the simulation results, which shows a good agreement. For the second one, the simulation results from three different models are compared, for both the local field distributions and the global impedances. It is shown that the EQS model can be used as an indicator to know at which frequency the Darwin model should be applied.

**Keywords:** electromagnetic; finite-element method; quasistatic models; resistive; capacitive; inductive effects

## 1. Introduction

The modeling of electrical devices is, at present, mainly based on the eddy current problem which allows a representative model of the behavior in the low-frequency domain to be obtained. With the advent of power electronics, devices are subject to high-frequency voltage and current stresses that must be taken into account, since this leads to the accelerated aging of insulators. Meanwhile, due to the development of the wide bandgap semiconductors used in the power converters, the frequency of voltage waveforms applied to the winding of electrical devices becomes increasingly high, especially in the case of pulse-width modulation (PWM). Consequently, resistive, capacitive, and even inductive effects should be simultaneously considered.

In the literature, the electrostatic, electrokinetic, and magnetostatic problems are widely used at low frequencies to capture the decoupled capacitive, resistive, and inductive effects, respectively [1,2]. When the frequency of the voltage supply applied to windings of electrical devices increases, the electric field in the regions close to the windings cannot be neglected, especially in the case of PWM [3]. In this case, special attention should be paid to the coupled resistive-capacitive effects. In practice, the electro-quasistatic (EQS) model [4] is commonly used for fields resulting from high voltage, withstand applications, or microelectronics, since this model is valid when the characteristic length of magnetic phenomena is low with respect to the wavelength [4,5]. However, if the skin effect appears

<footnote n="1">¹ Univ. Lille, Arts et Metiers Institute of Technology, Centrale Lille, Junia, ULR2697-L2EP, 59000 Lille, France; houssein.taha@univ-lille.fr (H.T.); thomas.henneron@univ-lille.fr (T.H.); yvonnick.le-menach@univ-lille.fr (Y.L.M.); florentin.salomez@univ-lille.fr (F.S.)</footnote>

<footnote n="2">² EDF R&D, ERMES, 7 Boulevard Gaspard Monge, 91120 Palaiseau, France; jean-pierre.ducreux@edf.fr</footnote>

<footnote>* Correspondence: zuqi.tang@univ-lille.fr</footnote>
on the distribution of the current density in the winding, the coupled resistive-inductive effects have to be taken into account, and so the magneto-quasistatic (MQS) [1] model is generally used.

Moreover, if all these coupled effects are characterized at the same time, the full Maxwell model [2,6] can be used, even it is more time-consuming and has an instability issue at low frequencies, as reported in [6,7]. On the other hand, if the modeling is only considered in the intermediate frequency range, the Darwin model can be a good choice [8–12], which includes all the capacitive, resistive, and inductive effects, but neglects the radiation one.

In the literature, the classification or the theoretical limitations of the different quasistatic models are mainly based on theoretical considerations, with several assumptions for the computational domain. For example, only the conductive domain is considered, not the case with multi-domains. The existing reference [13] showed a comparison of the results obtained with different quasistatic models in the time domain, when a simple axisymmetric test model, represented by a parallel plate capacitor, is used. The comparison results were illustrated for the electric field.

The aim of our work is to investigate the limits of the above-mentioned quasistatic models in the general case, particularly for the computational domain, including both conductive and non-conductive domains, which has not been well addressed in the literature. We are interested in the frequency domain due to the impedances' computation.

To achieve this objective, in this work, we first validated the Darwin model with an industrial example by comparing our numerical results and measurement results. To the best of the author's knowledge, the comparison of the Darwin model with real measurement results has not been reported in the literature. Secondly, with an academic example, we compared all the quasistatic models in the frequency domain, by regarding different numerical aspects, namely, the current density, the electric field, the magnetic field, and the impedance.

The paper is organized as follows: the full-wave Maxwell system is briefly introduced in Section 2. In Section 3, the EQS, MQS, and Darwin models are presented. In Section 4, two different examples, an industrial application and an academic electromagnetic device are provided for the numerical parts. Finally, the conclusion is given in Section 5.

## 2. Maxwell's System

Let us consider $\Omega \subset \mathbb{R}^3$, an open connected domain with boundary $\Gamma = \Gamma_H \cup \Gamma_B = \partial\Omega$ and $\Gamma_H \cap \Gamma_B = \emptyset$. $\Omega_c \subset \Omega$ presents the conductive sub-domain and $\Omega_{nc} \subset \Omega$ is the non-conductive sub-domain, as shown in Figure 1.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** The domain of study.</figcaption></figure>
The full Maxwell system, listed below, describes the magnetic and electric behaviors of the electromagnetic device

$$ \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}, \qquad (1) $$

$$ \nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}, \qquad (2) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (3) $$

$$ \nabla \cdot \mathbf{D} = 0 \qquad (4) $$

associated with the following boundary conditions

$$ \Gamma_E : \mathbf{E} \times \mathbf{n} = \mathbf{0}, \qquad (5) $$

$$ \Gamma_B : \mathbf{B} \cdot \mathbf{n} = 0, \qquad (6) $$

$$ \Gamma_H : \mathbf{H} \times \mathbf{n} = \mathbf{0}, \qquad (7) $$

$$ \Gamma_J : \mathbf{J} \cdot \mathbf{n} = 0 \qquad (8) $$

followed by the behavior laws of materials

$$ \mathbf{B} = \mu \mathbf{H}, \qquad (9) $$

$$ \mathbf{J} = \sigma \mathbf{E}, \qquad (10) $$

$$ \mathbf{D} = \varepsilon \mathbf{E} \qquad (11) $$

where E is the electric field (V/m), H is the magnetic field (A/m), B is the magnetic flux density (T), J is the current density (A/m²), D is the electric flux density (C/m²), μ is the magnetic permeability, σ and ε are the electric conductivity and permittivity, respectively.

### 3. Different Quasistatic Models

In low frequencies, by neglecting the radiation effects, several approximated models can be derived from Maxwell's equations based on the known Galilean limits [14,15]. There are three main different approaches: one is electric, called the EQS model, one is magnetic, called the MQS model, and one is electric and magnetic at the same time, called the Darwin model.

Let us denote v the velocity of the system with modulus |v| = L/T, where L and T denote the units of space and time, respectively. We compare this with the light celerity in the continuous medium c = 1/√εμ. In the literature, the condition |v| ≪ c is used to classify the quasistatic models, but this condition is not sufficient to distinguish between the EQS, MQS, and Darwin models.

In vacuum, if we consider this in the one-dimensional case, Equation (1) directly implies

$$ \frac{\partial E'}{\partial x} \sim \frac{\partial B'}{\partial t} \qquad (12) $$

then, we can obtain

$$ \frac{E'}{L} \sim \frac{B'}{T} \Rightarrow E' \sim |v|B' \qquad (13) $$

where X' is the order of magnitude of X, (X represents the electromagnetic fields E or B). The notation x ~ y means that the quantities x and y have the same order of magnitude.

Similarly, from Equation (2), we can obtain

$$ \frac{\partial H'}{\partial x} \sim \frac{\partial D'}{\partial t} \Rightarrow \frac{B'}{\mu L} \sim \frac{\varepsilon E'}{T} \Rightarrow B' \sim \frac{|v|}{c^2} E'. \qquad (14) $$
### 3.1. EQS Model

Suppose that assumption (14) holds and $E' \gg |v|B'$, the EQS model has to be adopted. This means that the term $j\omega\mathbf{B}$ defined in the Maxwell-Faraday's law, as given in Equation (1), becomes negligible, which implies that the induced current density is neglected. Applying the divergence operator to Equation (2), the charge conservation law in the frequency domain reads

$$ \nabla \cdot (\mathbf{J} + j\omega\mathbf{D}) = 0. \qquad (15) $$

Then, an electric scalar potential (ESP) $\varphi$ can be introduced, such as

$$ \mathbf{E} = -\nabla\varphi \text{ in } \Omega, \qquad (16) $$

$$ \varphi = c_1 \text{ on } \Gamma_E \qquad (17) $$

with $c_1$ as a constant. By combining the expression of $\mathbf{E}$ with (10), (11), and (15), the potential formulation for the EQS model reads

$$ \nabla \cdot (\sigma \nabla \varphi + j\omega\varepsilon \nabla \varphi) = 0 \text{ in } \Omega. \qquad (18) $$

To apply a voltage excitation on the terminals of the winding corresponding to the sub-domain $\Omega_c$, Dirichlet boundary conditions are imposed on $\varphi$. By applying the finite element (FE) method, the ESP $\varphi$ is discretized using the nodal elements, which introduces a symmetric matrix system.

### 3.2. MQS Model

Suppose that assumption (13) holds and $B' \gg |v|/c^2E'$: the MQS model has to be adopted. This means that the term $j\omega\mathbf{D}$ defined in the Maxwell-Ampere's law, as given in Equation (2), becomes negligible, which implies that the displacement current is neglected in the dielectrics. The MQS model is derived as follows

$$ \nabla \times \mathbf{E} = -j\omega\mathbf{B}, \qquad (19) $$

$$ \nabla \times \mathbf{H} = \mathbf{J}, \qquad (20) $$

$$ \nabla \cdot \mathbf{B} = 0. \qquad (21) $$

From (21), the magnetic vector potential (MVP) $\mathbf{A}$ can be introduced, such as

$$ \mathbf{B} = \nabla \times \mathbf{A}, \qquad (22) $$

$$ \mathbf{A} \times \mathbf{n} = 0 \text{ on } \Gamma_B \text{ and } \Gamma_E \qquad (23) $$

and the ESP $\varphi$ is only defined in the conductive sub-domain, as follows

$$ \mathbf{E} = -j\omega\mathbf{A} - \nabla\varphi, \qquad (24) $$

$$ \varphi = c_1 \text{ on } \Gamma_E. \qquad (25) $$

Combining the expressions of $\mathbf{E}$ and $\mathbf{B}$ with (9), (10), and (20), the $\mathbf{A}-\varphi$ potential formulation for the MQS model reads

$$ \nabla \times (\nu \nabla \times \mathbf{A}) + \sigma(j\omega\mathbf{A} + \nabla\varphi) = \mathbf{0} \text{ in } \Omega, \qquad (26) $$

$$ \nabla \cdot (\sigma(j\omega\mathbf{A} + \nabla\varphi)) = 0 \text{ in } \Omega_c \qquad (27) $$

where $\nu$ is the magnetic reluctivity (the inverse of the magnetic permeability $\mu$). Applying the FE method, the MVP $\mathbf{A}$ is discretized with the edge elements and the ESP $\varphi$ with the nodal elements, which introduce a symmetric matrix. It should be mentioned here that the resultant matrix is not gauged. The uniqueness of $\mathbf{A}$ can be ensured by adding the Coulomb gauge [15] or the tree gauge [16].
### 3.3. Darwin Model

By the Helmholtz decomposition, the electric field **E** can be split into two parts: an irrotational part **E**<sub>irr</sub>, which is curl free, and a solenoidal part **E**<sub>sol</sub>, which is divergence free [17]

$$ \mathbf{E} = \mathbf{E}_{\text{irr}} + \mathbf{E}_{\text{sol}}, \quad \nabla \times \mathbf{E}_{\text{irr}} = \mathbf{0}, \quad \nabla \cdot \mathbf{E}_{\text{sol}} = 0. \qquad (28) $$

We rewrite the Maxwell equations using this decomposition in vacuum

$$ \nabla \times \mathbf{E}_{\text{sol}} = -j\omega\mathbf{B}, \qquad (29) $$

$$ \nabla \times \mathbf{H} = j\omega\varepsilon\mathbf{E}_{\text{irr}} + j\omega\varepsilon\mathbf{E}_{\text{sol}}, \qquad (30) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (31) $$

$$ \varepsilon \nabla \cdot \mathbf{E}_{\text{irr}} = 0. \qquad (32) $$

Based on the above system and on the Galilean limits, there are two orderings [18]

$$ \frac{E'_{\text{sol}}}{B'} \sim |v| \quad \text{and} \quad \frac{E'_{\text{irr}}}{B'} \sim \frac{1}{|v|}. \qquad (33) $$

This mixed ordering eliminates radiation effects; then, the Darwin model reads

$$ \nabla \times \mathbf{E} = -j\omega\mathbf{B}, \qquad (34) $$

$$ \nabla \times \mathbf{H} = \mathbf{J} + j\omega\varepsilon\mathbf{E}_{\text{irr}}, \qquad (35) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (36) $$

$$ \nabla \cdot \mathbf{D} = 0. \qquad (37) $$

It can be observed that the Darwin model neglects the rotational part of the displacement current densities with $j\omega\varepsilon\mathbf{E}_{\text{rot}} = 0$. From (28), $\mathbf{E}_{\text{irr}}$ can be represented by the gradient of the scalar electric $\varphi$ while $\mathbf{E}_{\text{sol}}$ can be represented by the time derivative of the magnetic vector

$$ \mathbf{E} = \mathbf{E}_{\text{irr}} + \mathbf{E}_{\text{sol}} = -\nabla\varphi - j\omega\mathbf{A}. \qquad (38) $$

Combining (38) with (9), (10), and (35), the $\mathbf{A}-\varphi$ potential formulation for the Darwin model [9] reads

$$ \begin{aligned} & \nabla \times (\nu \nabla \times \mathbf{A}) + \sigma(j\omega\mathbf{A} + \nabla\varphi) + j\omega(\varepsilon\nabla\varphi) = \mathbf{0}, \\ & \nabla \cdot (\sigma(j\omega\mathbf{A} + \nabla\varphi) + j\omega(\varepsilon\nabla\varphi)) = 0. \end{aligned} \qquad (39) $$

It should be mentioned here that the ESP $\varphi$ is defined in the whole domain and the boundary conditions are the same as for the MQS model.

However, differing from the case of MQS model, the resultant FE matrix is non-symmetric. Moreover, the matrix is ill-conditioned due to the high contrast in the magnitudes of material coefficients [9], which makes the convergence rate very slow, with the iterative solvers. A Coulomb type gauge as $\nabla \cdot (\varepsilon\mathbf{A}) = 0$ can be imposed to make the matrix system symmetric, and direct solvers are used to ensure the stability of the convergence [11]. However, it is difficult to use the direct solvers for industrial applications due to the memory limitations, since the number of degrees of freedom (DoF) is relatively large in practice. Thus, in this work, the iterative solver is always preferred. Without the preconditioner, the convergence with iterative solvers is sometimes prone to errors, and causes a delay in convergence and failure to find an acceptable solution.

In our study, the biconjugate gradient stabilized (BiCGSTAB) solver [19] is implemented with a Split-Jacobi preconditioner to solve the unsymmetric system $\mathbf{A}\mathbf{x} = \mathbf{b}$ given in (39). The deriving preconditioner splits the Jacobi preconditioner $\mathcal{M}$ ($\mathcal{M}$ is the diagonal of $\mathbf{A}$) as $\mathcal{M} = \mathcal{S}^T\mathcal{S}$, where $\mathcal{S}_{i,j} = \sqrt{\mathcal{M}_{i,j}}$. Then, the final system to solve is given as follows

$$ \mathcal{S}^{-1}\mathbf{A}(\mathcal{S}^T)^{-1}(\mathcal{S}^T\mathbf{x}) = \mathcal{S}^{-1}\mathbf{b}. \qquad (40) $$
## 4. Numerical Validation

All the formulations presented in Section 3 were implemented in our academic software *code_Carmel* (https://code-carmel.univ-lille.fr/, accessed on 20 November 2021), which provided the numerical simulation results. In the following, two different examples were carried out. For the first example, we validated the simulation results obtained with the Darwin model by considering an industrial transformer. The experimental measurement data, in a common mode test, were provided until 100 MHz. However, based on the measurement results, it should be mentioned here that the Darwin model is not valid when the frequency is superior to 10 MHz, since the radiation effects become non-neglected. As the frequency increases, to handle the skin effect in the periphery of the windings, more elements are needed in the mesh, which introduces a scientific challenge at the numerical aspect. For our computation, the total DoFs reached 40 million elements. The second one, which represents an inductance example, is studied to investigate the influence of the frequency on the electromagnetic field distributions and compare the impedance curves computed by each model, namely, EQS, MQS, and Darwin models.

### 4.1. Industrial Application

For the first example, the geometry consists of two windings wound in the same direction around a toroidal core, as shown in Figures 2–4. The conductor used for the windings was made of copper, presented in red in Figure 4: it was circular in shape and each conductor was composed of 10 turns with a section of 0.5027 mm². The enamel component covered the winding with a thickness of 40 µm. The toroid used for this experiment was a TDK-Epcos toroid reference B64290L0618X830 (https://product.tdk.com/en/search/ferrite/ferrite/ferrite-core/info?part_no=B64290L0618X830, accessed on 20 November 2021), of material N30 presented in gray in Figure 3, where the material in blue is the embedding with a thickness of 0.25 mm. The conductivity of the conductor is equal to 58 MS/m, while the enamel was emphasized so that the relative permittivity was set as 5.4.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Toroidal core with two windings-side.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** 3-D model of the transformer.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Detail of winding.</figcaption></figure>

### 4.1.1. Computational Configurations

An electric potential difference $φ = V_s$ was imposed between the terminals of the conductor of the windings marked in red, as shown in Figure 4. A single-phase common mode test was considered in this example, which means that $φ = V_s$ was imposed on one terminal of the conductor and $φ = -V_s$ on the other, as shown in Figure 3. The frequency interval was [0:10⁶] Hz. To handle the skin effect, the considered mesh was well-refined and suitable for all frequencies up to 1 MHz, which featured 29,832,477 tetrahedrons including 5,257,323 nodes and 35,093,396 edges. For each frequency, the computing time using the Darwin model took about 9 days for 40,338,477 DoFs using the BiCGSTAB solver. Notably, for a system with a large DoFs, the BiCGSTAB is always preferred to direct solvers because the latter may fail due to the relative need for memory, which is much larger than that required for the iterative solvers. The mesh of the winding and the toroidal core are presented in Figure 5.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Considered mesh.</figcaption></figure>

### 4.1.2. Electric Potential and Electric Field

The electric scalar potential along the magnetic core is presented in Figure 6a, which shows a decreasing linear variation in its magnitude. In addition, the electric field distribution in a sectional view is also given in the magnetic core, as shown in Figure 6b. It can be observed that the electric field is located in the periphery of the core due to the skin effect that occurs, and reaches a maximum value of 12 V/m.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Distribution of the electric scalar potential and the electric field in the magnetic core for $f = 500$ kHz.</figcaption></figure>

### 4.1.3. Conduction and Displacement Current Densities

When the frequency increases, the total current is divided into two parts: the conduction current in the conductive domains and the displacement current located in the non-conductive domains. This means that the displacement current intervenes in the modeling of the electric field, in addition to the eddy current.

On one hand, the conduction current density is presented in Figure 7a in the windings for a frequency of 500 kHz. High values are observed inside each of the two conductors. Due to the skin effect, the current is concentrated in the small layer at the boundary of the conductor. Moreover, the peaks were higher near the middle section, mainly due to the magnetic coupling known as the proximity effect.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Distribution of the conduction current density and the displacement current density in the windings and the magnetic core, respectively.</figcaption></figure>

On the other hand, as shown in Figure 7b, the displacement current density is presented in a sectional view of the magnetic core for a frequency of 500 kHz. The distribution
shows that the significant quantity of current is located between the turns, as the magnitude tends to be zero in the contact area with the windings, due to the capacitive effects.

#### 4.1.4. Measurement and Simulation Comparison

In Figure 8, the evolution of the impedance as a function of the frequency obtained from Darwin model, as well as the phase, are presented.

First, the impedance Z corresponds to the DC resistance of the winding when the frequency tends towards zero. Indeed, we have $Z = R_{DC} = 6.8$ mΩ.

Second, for $f < 500$ kHz, the results obtained with the Darwin are in good agreement with the measurement results, which present the resistive-inductive phenomena, such as the MQS regime.

Third, the resonant frequency is around 1 MHz; the simulation result of the Darwin model shows a close result at the chosen point of the frequency (7th point), just before the resonant frequency. It can be found that this point is no longer in the MQS regime, since the impedance at this point is not linear to the previous points. A similar behavior can be observed for the phase as a function of the frequency.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Modulus of impedance and the phase computed with respect to frequencies.</figcaption></figure>

### 4.2. Modeling of the Academic Electromagnetic Device

For the second example, we are interested in determining the operating frequency domain of different models. Some validation work was carried out with the static or/and MQS and EQS models, as shown in [20,21]. The objective was to study the influence of the frequency on the distribution of electromagnetics fields, as well as on the impedance obtained using the models mentioned in Section 3.

As shown in Figure 9, the electromagnetic device composed of a winding (shown in red color) coiled around the central column of the magnetic core (shown in blue color) was considered. We have a conductor of section 11.36 mm² which takes a rectangular shape. It is composed of 11 turns and has a length of 2.3 m. The insulation covered the winding with a thickness of 0.25 mm. The magnetic core was composed of N30-ferrite. This inductance is shown in 2-D cutting plane in Figure 10.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Geometry of the electromagnetic device.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Schematic representation of the electromagnetic device.</figcaption></figure>

The conductivity of the copper coil is set as 58 MS/m. The insulation covered the winding, and it is emphasized that the relative permittivity was set as 9.4 when the magnetic core is composed of N30-ferrite, characterized by a relative permeability of 4300, an electric permittivity of 80 and an electric conductivity fixed to zero. The electric potential on the below face of the magnetic core was taken to be zero, which corresponds to connect to the ground.

### 4.2.1. Computational Configurations

A sinusoidal voltage was applied between the terminals of the winding marked in red, as shown in Figure 9. The frequency interval was [0:10⁵] Hz. To consider the skin effects in the winding for high frequencies, the mesh used for the model featured 14,443,563 tetrahedrons, including 2,466,531 nodes and 16,911,563 edges. The mesh of the winding and the magnetic core are presented, respectively, in Figure 11a,b. For the EQS model, the number of DoFs is 2,466,002 and the computational time for one frequency takes about two hours, while it takes about 75 hours for 19,316,839 DoFs in the MQS model. The computing time using Darwin formulation takes about 86 hours for 19,373,155 DoFs. It should be recalled here that the matrix system associated with MQS is symmetric, while this is not the case for the Darwin model. In the general case, the computational cost for the MQS model is less expensive than Darwin's due to the ill-conditioned matrix of the Darwin model, and different preconditioners can be applied for MQS and Darwin cases.
Moreover, the computational cost for the EQS model is relatively low in comparison with the two other models.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Considered mesh.</figcaption></figure>

Additionally, It should be mentioned that the condition number of the matrix depends on the frequency for both MQS and Darwin models. Consequently, the number of iterations associated with the iterative solver, which is required to achieve the solution's convergence, is not the same. For example, for *f* = 10 kHz and with a fixed stopping criterion, the number of iterations associated with BiCGSTAB is 14 k using the MQS model, while it takes about 17 k for the Darwin model, as shown in Table 1.

<table><caption>**Table 1.** Comparison between the different models for a frequency *f* = 10 kHz.</caption><thead><tr><td>Formulations</td><td>DoFs</td><td>Iterations</td><td>CPU Time</td></tr></thead><tbody><tr><td>EQS</td><td>2,466,002</td><td>18 k</td><td>≈2.5 h</td></tr><tr><td>MQS</td><td>19,316,839</td><td>14 k</td><td>≈75 h</td></tr><tr><td>Darwin</td><td>19,373,155</td><td>17 k</td><td>≈86 h</td></tr></tbody></table>

### 4.2.2. Distribution of Electric and Magnetic Fields

At low frequencies, when the frequency tends towards zero, the distribution of the current density in the winding obtained with all models is similar, as observed in Figure 12a, since this is equivalent to solve an electrokinetic problem. In this case, the distribution of the current density is only influenced by the geometry of the winding. When the frequency supplying the winding increases, coupled effects may be observed. To represent these effects in the winding, the distribution of the current density *J* is presented in Figure 12b-d. At *f* = 10 kHz, with the EQS model, it is clear that *J* will never be homogeneous along the winding; the majority of currents flow through the capacitor in terms of the displacement current. On the other hand, the current density with the MQS and Darwin models is restricted to a small layer at the boundary of the winding, and starts to flow in the outskirts of the conductor, as observed by the skin effect.

Besides, due to the capacitive effects, the current density distribution obtained with the Darwin model is not the same as that of MQS. The irrotational current density is much larger in both cases, compared to the solenoidal one, such that the capacitive coupling cannot be distinguished from Figure 12d. Moreover, the distribution of the electric field obtained from the Darwin model is presented in Figure 13 in a 2-D cutting plane for *f* = 10 kHz. A high field strength appears in the dielectric regions, particularly between the gaps of the winding.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Distribution of the current density *J* (A/m²) in the winding.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Distribution of the electric field *E* (V/m) between the gaps of the winding at *f* = 10 kHz.</figcaption></figure>

In addition, when capacitive effects are neglected for low and middle frequencies, both MQS and Darwin problems give the same distribution of the magnetic flux density. Then, when the frequency tends towards zero, the distribution of the magnetic flux density
obtained from MQS and Darwin models is equivalent to solve a magnetostatic problem, as shown in Figure 14.

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Distribution of the magnetic flux density **B** (T) in the magnetic core at *f* = 1 Hz.</figcaption></figure>

4.2.3. Evolution of the Impedance Versus the Frequency

The evolution of the impedance modulus and its phase, according to the frequency obtained from EQS, MQS, and Darwin models, are presented in Figure 15. At low frequencies, the impedance Z corresponds to the DC resistance of the winding. Then, for all models, we obtain Z = RDC = 3 mΩ when the frequency tends towards zero. When the frequency increases, the skin effect appears in the winding; then, the value of the resistance of the winding increases. This effect is not taken into account with the EQS model. For f ∈ [0:10³] Hz, the MQS and Darwin models give a similar evolution of the impedance. Indeed, the capacitive effect is negligible compared with the inductive and resistive effects. For f > 1 kHz, the influence of the capacitive effects appears on the evolution of the impedance for both the EQS and Darwin models. A similar behavior can be observed for both models in Figure 15.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Modulus of impedance and the phase computed with respect to frequencies.</figcaption></figure>

**5. Conclusions**

In this work, the measurement results for an industrial example have been provided to validate the simulation results obtained with the Darwin model, which shows a good agreement in the range of intermediate frequencies, particularly around the resonant frequency. Additionally, for a complex and multi-domain application, the different quasistatic models were compared by computing their impedances, as well as the different electromagnetic fields, to investigate their limits.

At low frequencies, when the impedance exhibits a resistive behavior, the EQS model is sufficient. Moreover, when the inductive effects are non-neglected, the MQS and Darwin models can be invoked, but the MQS model is preferred, since its computational cost is lower than the Darwin model. However, in the intermediate frequencies, particularly around the resonant frequency, the Darwin model should be adopted, since the MQS model cannot handle the coupled capacitive-inductive effects.

Furthermore, the EQS model can be used as an indicator to learn the frequency at which the capacitive effects become significant, where the Darwin model should be applied.
**Author Contributions:** All authors are involved to some extent in the different parts of the paper: Conceptualization, H.T., Z.T., T.H., Y.L.M. and J.-P.D.; Data curation, H.T.; Formal analysis, H.T., Z.T., T.H. and Y.L.M.; Funding acquisition, Y.L.M. and J.-P.D.; Investigation, H.T., Z.T. and T.H.; Methodology, Z.T., T.H. and Y.L.M.; Resources, F.S. and J.-P.D.; Software, H.T.; Supervision, Z.T., T.H., Y.L.M. and J.-P.D.; Validation, Z.T., T.H. and Y.L.M.; Visualization, H.T., Z.T. and T.H.; Writing—original draft, H.T.; Writing—review & editing, Z.T. and T.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been financially supported by the European Regional Development Fund (ERDF) NP0017752 and LAMEL laboratory (joint laboratory between the University of Lille, ENSAM, and the EDF company).

**Data Availability Statement:** The data presented in this study are available on request from the
corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interests.

## References

1. Haus, H.; Melcher, J. *Electromagnetic Fields and Energy*; Prentice Hall: Hoboken, NJ, USA, 1989.

2. Badics, Z.; Pávó, J. Full Wave Potential Formulation with Low-Frequency Stability Including Ohmic Losses. *IEEE Trans. Magn.* **2015**, *51*, 1–4. [CrossRef]

3. Maki-Ontto, P.; Luomi, J. Induction motor model for the analysis of capacitive and induced shaft voltages. In Proceedings of the IEEE International Conference on Electric Machines and Drives, San Antonio, TX, USA, 15 May 2005; pp. 1653–1660. [CrossRef]

4. Steinmetz, T.; Helias, M.; Wimmer, G.; Fichte, L.O.; Clemens, M. Electro-quasistatic field simulations based on a discrete electromagnetism formulation. *IEEE Trans. Magn.* **2006**, *42*, 755–758. [CrossRef]

5. Weida, D.; Richter, C.; Clemens, M. Design of ZnO microvaristor material stress-cone for cable accessories. *IEEE Trans. Dielectr. Electr. Insul.* **2011**, *18*, 1262–1267. [CrossRef]

6. Hiptmair, R.; Krämer, F.; Ostrowski, J. A Robust Maxwell Formulation for All Frequencies. *IEEE Trans. Magn.* **2008**, *44*, 682–685. [CrossRef]

7. Bebendorf, M.; Krämer, F. Hierarchical Matrix Preconditioning for Low-Frequency-Full-Maxwell Simulations. *Proc. IEEE 2013*, **101**, 423–433. [CrossRef]

8. Larsson, J. Electromagnetics from a quasistatic perspective. *Am. J. Phys.* **2006**, *75*, 230–239. [CrossRef]

9. Koch, S.; Schneider, H.; Weiland, T. A Low-Frequency Approximation to the Maxwell Equations Simultaneously Considering Inductive and Capacitive Phenomena. *IEEE Trans. Magn.* **2012**, *48*, 511–514. [CrossRef]

10. Eller, M.; Reitzinger, S.; Schöps, S.; Zaglmayr, S. A Symmetric Low-Frequency Stable Broadband Maxwell Formulation for Industrial Applications. *SIAM J. Sci. Comput.* **2017**, *39*, B703–B731. [CrossRef]

11. Zhao, Y.; Tang, Z. A Novel Gauged Potential Formulation for 3-D Electromagnetic Field Analysis Including Both Inductive and Capacitive Effects. *IEEE Trans. Magn.* **2019**, *55*, 1. [CrossRef]

12. Clemens, M.; Kähne, B.; Schöps, S. A Darwin Time Domain Scheme for the Simulation of Transient Quasistatic Electromagnetic Fields Including Resistive, Capacitive and Inductive Effects. In Proceedings of the 2019 Kleinheubach Conference, Miltenberg, Germany, 23–25 September 2019; pp. 1–4.

13. Koch, S.; Weiland, T. Different types of quasistationary formulations for time domain simulations. *Radio Sci.* **2011**, *46*, 1–7. [CrossRef]

14. Rapetti, F.; Rousseaux, G. Implications of Galilean electromagnetism in numerical modeling. In Proceedings of the IET 8th International Conference on Computation in Electromagnetics (CEM 2011), Wroclaw, Poland, 11–14 April 2011; pp. 1–2. [CrossRef]

15. Bellac, M.; Levy-Leblond, J.M. Galilean electromagnetism. *Il Nuovo Cim. B* **1973**, *14*, 217–234. [CrossRef]

16. Albanese, R.; Rubinacci, G. Magnetostatic field computation in terms of two-component vector potentials. *Int. J. Numer. Methods Eng.* **1990**, *29*, 515–532. [CrossRef]

17. Cortes Garcia, I.; Schöps, S.; De Gersem, H.; Baumanns, S. Systems of Differential Algebraic Equations in Computational Electromagnetics. In *Differential-Algebraic Equations Forum*; Springer: Cham, Switzerland, 2018; pp. 123–169. [CrossRef]

18. Kruger, S.E. The Three Quasistatic Limits of the Maxwell Equations. *arXiv* **2019**, arXiv:1909.11264

19. Van der Vorst, H.A. Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems. *SIAM J. Sci. Comput.* **1992**, *13*, 631–644. [CrossRef]

20. Steinmetz, T.; Kurz, S.; Clemens, M. Domains of Validity of Quasistatic and Quasistationary Field Approximations. In Proceedings of the VXV International Symposium on Theoretical Engineering, Lubeck, Germany, 22–24 June 2009; pp. 1237–1247.

21. Mazauric, V.G.; Rondot, L.; Wendling, P.F. Enhancing Quasi-Static Modeling: A Claim for Electric Field Computation. *IEEE Trans. Magn.* **2013**, *49*, 1629–1632. [CrossRef]