*4.4. Final System of Equations*

Equations (6), (20) and (25) are solved by using a numerical method. The most suitable way is to mesh the averaged surface Γ into *n* triangular elements and *p* nodes. Let us assume that the tangential component of the eddy current and the magnetization in each element are uniform.

*Appl. Sci.* **2020**, *10*, 4284

The Galerkine projection method is then applied to Equation (6), and we get the following matrix system [11,15,16]:

$$
\begin{bmatrix} \mathbf{A} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{M} \end{bmatrix} + \begin{bmatrix} \mathbf{B} \end{bmatrix} \cdot \begin{bmatrix} \Delta \Phi \end{bmatrix} = 0,\tag{26}
$$

where:

[**M**] is a complex vector of dimension 3*n*; [Δ**Φ**] is a complex vector of dimension *p*;

[**A**] is a (*p* × *3n*) matrix expressed as follows:

$$\mathbf{A}(i,k) = 2j\omega \frac{\mu\_0 \mu\_r}{\mu\_r - 1} \int\_{\Gamma} w\_i \cdot \mathbf{n}\_k \mathbf{d}\Gamma = 0,\tag{27}$$

[**B**] is a (*p* × *p*) sparse matrix and can be written as:

$$\mathbf{B}(i,k) = (\alpha + \beta) \int\_{\Gamma} \mathbf{grad}\_{\delta} w\_{i} \cdot \mathbf{grad}\_{\delta} w\_{k} d\Gamma. \tag{28}$$

A matrix system representing Equation (20) can also be determined thanks to a point matching approach at the element centroids. This linear matrix system is expressed as [10,11,16]:

$$
\left(\frac{[\mathbf{I}\_d]}{\mu\_I - 1} + [\mathbf{F}]\right) \cdot [\mathbf{M}] - [\mathbf{C}] \cdot [\Delta \Phi] - [\mathbf{D}] \cdot [\mathbf{I}] = \mathbf{0},\tag{29}
$$

where [**I***d*] represents the identity matrix; [**F**] is a (3*n* × 3*n*) dense matrix; [**C**] is a (3*n* × *p*) also dense matrix; [**D**] is a (3*n* × *m*) matrix and [**I**] is a (*m* × 1) vector. The coefficients of these matrices are expressed as follows:

$$F(i, \text{ } j) = \left[ \left( \frac{\overline{\mathbf{G}}}{4\pi} \mathbf{grad} \, \int\_{\Gamma\_j} \frac{\mathbf{u}\_j \cdot \mathbf{r}\_i}{\mathbf{r}\_i^3} d\Gamma \right), \mathbf{u}\_i \right] \tag{30}$$

where [,] is defined as the scalar product operator and **u***<sup>i</sup>* is the vector basis of the element *i* and **r**<sup>i</sup> is the vector between the integration point on Γ<sup>j</sup> to the centroid of the element *i*.

$$\mathbf{C}(i,j) = \frac{1}{4\pi} \int\_{\Gamma} \frac{\mathbf{n}\_{\rangle} \times \mathbf{grad} w\_{i} \times \mathbf{r}\_{i}}{\mathbf{r}\_{i}^{3}} d\Gamma. \tag{31}$$

$$\mathbf{D}(i,j) = \frac{1}{4\pi} \frac{1}{\mathbf{S}\_j} \int\_{\Omega \mathbf{c}\_j} \frac{\mathbf{l}\_j \times \mathbf{r}\_i}{\mathbf{r}\_i^3} d\Omega. \tag{32}$$

By integrating Equation (25) on each conductor, we have the link between the partial voltages of the conductors to the currents flowing in them. Thus, the voltage appearing on the conductor *k* is given by:

$$\mathbf{U}\_{k} = \mathbf{R}\_{k}\mathbf{I}\_{k} + j\omega \sum\_{i=1}^{m} m\_{ik}\mathbf{I}\_{i} + j\omega \sum\_{i=1}^{n} m\_{ik}^{f}\mathbf{M}\_{\mathbf{a}i} + j\omega \sum\_{i=1}^{p} m\_{ik}^{q} \cdot \phi\_{i\nu} \tag{33}$$

where R*<sup>k</sup>* is the resistance of the *kth* conductor; m*ik* is the mutual inductance between the two conductors *k*, *i*; *m<sup>f</sup> ik* is defined as the mutual inductance between the magnetization **M**a*<sup>i</sup>* of the shell element and the conductor *k*; *mq ik* is defined as the mutual inductance between the scalar magnetic potential discontinuity Δφ of the shell node and the conductor *k*. These mutual inductances are expressed as follows:

$$m\_{i\bar{k}} = \frac{\mu\_0}{4\pi} \frac{1}{S\_i S\_k} \int\_{\Omega c\_k} \int\_{\Omega c\_i} \frac{\mathbf{l}\_i \mathbf{l}\_k}{\mathbf{r}} d\Omega\_{\mathrm{ci}} d\Omega\_{\mathrm{ck}} \tag{34}$$

$$m\_{ik}^f = \frac{\mu\_0 \mathbf{G}}{4\pi} \frac{1}{S\_k} \int\_{\Omega \mathbf{x}\_k} \int\_{\Gamma\_i} \frac{\mathbf{u}\_i \times \mathbf{r}}{\mathbf{r}^3} d\Gamma \mathbf{l}\_k d\Omega\_{ck\prime} \tag{35}$$

$$m\_{\rm ik}^{q} = \frac{\mu\_{0}}{4\pi} \frac{1}{S\_{\rm k}} \int\_{\Omega c\_{\rm k}} \int\_{\Gamma\_{i}} \frac{\mathbf{n}\_{i} \times \mathbf{grad} w\_{i}}{\mathbf{r}} d\Gamma \mathbf{l}\_{k} d\Omega\_{ck\prime} \tag{36}$$

where **u***<sup>i</sup>* is a basis vector of magnetizations; **n***<sup>i</sup>* is the normal vector of the element Γ*<sup>i</sup>* and **l***<sup>k</sup>* is the vector unit of the current direction in the conductor *k.* Let us note that the mutual inductance terms <sup>m</sup>*ik* can be calculated by the well-known semi-analytic formulas in [7–9,22,23]. Whereas, the terms *<sup>m</sup><sup>f</sup> ik*, *mq ik* are computed with a standard numerical integration.

Writing Equation (33) for all conductors, we get a matrix system known as the impedance matrix system:

$$\mathbf{[U]} = [\mathbf{Z}] \cdot [\mathbf{I}] + \left[\mathbf{L}^f\right] \cdot [\mathbf{M}] + [\mathbf{L}^g] \cdot [\boldsymbol{\Delta} \boldsymbol{\Phi}],\tag{37}$$

where [**Z**] is a (*<sup>m</sup>* <sup>×</sup> *<sup>m</sup>*) matrix; [**U**] is a (*<sup>m</sup>* <sup>×</sup> 1) vector; <sup>1</sup> **L***f* 2 is a (*<sup>m</sup>* <sup>×</sup> <sup>3</sup>*n*) matrix and [**L***q*] is a (*<sup>m</sup>* <sup>×</sup> *<sup>p</sup>*) matrix.

Finally, the algebraic linear systems (26), (29) and (37) is considered as *p* + *3n* + *m* complex unknowns:

$$
\begin{bmatrix}
\mathbf{A} & \mathbf{B} & \mathbf{0} \\
\mathbf{M}\mathbf{0}\mathbf{M} & -\mathbf{C} & -\mathbf{D} \\
\mathbf{L}^f & \mathbf{L}^g & \mathbf{Z}
\end{bmatrix} \times \begin{bmatrix}
\mathbf{M} \\
\Delta\Phi \\
\mathbf{I}
\end{bmatrix} = \begin{bmatrix}
0 \\
0 \\
\mathbf{U}
\end{bmatrix} \tag{38}
$$

where **MoM** = [**I***d*] <sup>μ</sup>*r*−<sup>1</sup> <sup>+</sup> [**F**] is usually called the magnetic moment method matrix.

Let us note that Equation (29) has a disadvantage when the permeability of the material is high. Indeed, in such cases, the matrix term [**I***d*]/(μ*<sup>r</sup>* − 1) becomes very small in comparison with the [**F]** matrix. This can lead to the singularity of matrix [**MoM**]. It should be noted that for the linear case, the simpler formulations can be developed. In such cases, Laplace equation is applied in the whole electromagnetic thin shell regions. Formulation (30) can be determined by a surface integral equation generated by charge distributions (Coulombian approach) or current distributions (Amperian approach) located on the surface area bounding the volume Ω of the thin shell region. In such configuration, only this surface area must be discretized. Consequently, the degrees of freedom of the obtained matrix system is lower. Moreover, Newton-Raphson algorithm can be easily used for the non-linear case [25].

Besides, in order to solve the equation system (38) with a reduced number of degrees of freedom, a Kirchhoff's mesh rule is introduced. Consequently, the current and the voltage values become associated to each of the independent circuit mesh. The coupling formulation has been implemented for 3D geometry.

#### **5. Numerical Examples**

In this section, we consider two numerical examples. To valid our formulation, three numerical methods are compared. The first one is a 3D shell element FEM coupled with the circuit equations [5]. In this modelling, a special care is proposed for the mesh around the conductor to ensure accurate results. The second one is a 2D FEM formulation. This approach has already shown its good precision with a few numbers of elements and is considered as the reference. The last one is our coupling integral formulation.
