# Numerical Simulation-Based Investigation of the Limits of Different Quasistatic Models
## 3 Different Quasistatic Models
### Darwin Model  
By the Helmholtz decomposition, the electric field \(\mathbf{E}\) can be split into two parts: an irrotational part \(\mathbf{E}_{\text{irr}}\), which is curl free, and a solenoidal part \(\mathbf{E}_{\text{sol}}\), which is divergence free [17]  
\[\mathbf{E}=\mathbf{E}_{\text{irr}}+\mathbf{E}_{\text{sol}},\ \
abla\times\mathbf{E}_{\text{irr}}= \mathbf{0},\ \
abla\cdot\mathbf{E}_{\text{sol}}=0. \tag{28}\]  
We rewrite the Maxwell equations using this decomposition in vacuum  
\[\
abla\times\mathbf{E}_{\text{sol}} =-j\omega\mathbf{B}, \tag{29}\] \[\
abla\times\mathbf{H} =j\omega\epsilon\mathbf{E}_{\text{irr}}+j\omega\epsilon\mathbf{E}_ {\text{sol}},\] (30) \[\
abla\cdot\mathbf{B} =0,\] (31) \[\varepsilon\
abla\cdot\mathbf{E}_{\text{irr}} =0. \tag{32}\]  
Based on the above system and on the Galilean limits, there are two orderings [18]  
\[\frac{\text{E}^{\prime}_{\text{sol}}}{B^{\prime}}\sim|v|\quad\text{and}\quad \frac{\text{E}^{\prime}_{\text{irr}}}{B^{\prime}}\sim\frac{1}{|v|}. \tag{33}\]  
This mixed ordering eliminates radiation effects; then, the Darwin model reads  
\[\
abla\times\mathbf{E} =-j\omega\mathbf{B}, \tag{34}\] \[\
abla\times\mathbf{H} =\mathbf{J}+j\omega\epsilon\mathbf{E}_{\text{irr}},\] (35) \[\
abla\cdot\mathbf{B} =0,\] (36) \[\
abla\cdot\mathbf{D} =0. \tag{37}\]  
It can be observed that the [PERSON] model neglects the rotational part of the displacement current densities with \(j\omega\epsilon\mathbf{E}_{\text{rot}}=\mathbf{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}}=-\
abla\varphi-j \omega\mathbf{A}. \tag{38}\]  
Combining (38) with (9), (10), and (35), the \(\mathbf{A}\)-\(\varphi\) potential formulation for the [PERSON] model [9] reads  
\[\
abla\times(\
u\
abla\times\mathbf{A})+\sigma(j\omega\mathbf{A} +\
abla\varphi)+j\omega(\varepsilon\
abla\varphi) =\mathbf{0}, \tag{39}\] \[\
abla\cdot(\sigma(j\omega\mathbf{A}+\
abla\varphi)+j\omega( \varepsilon\
abla\varphi)) =0.\]  
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 \(\
abla\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 \(\mathbb{A}\mathbf{x}=\mathbf{b}\) given in (39). The deriving preconditioner splits the Jacobi preconditioner \(\mathbb{M}\) (\(\mathbb{M}\) is the diagonal of \(\mathbb{A}\)) as \(\mathbb{M}=\mathbb{S}^{T}\mathbb{S}\), where \(\mathbb{S}_{i,j}=\sqrt{\mathbb{M}_{i,j}}\). Then, the final system to solve is given as follows  
\[\mathbb{S}^{-1}\mathbb{A}(\mathbb{S}^{T})^{-1}(\mathbb{S}^{T}\mathbf{x})= \mathbb{S}^{-1}\mathbf{b}. \tag{40}\]