**1. Introduction**

All of the electromagnetic phenomena occurring in the electrical systems are described by Maxwell's equations, together with the constitutive material laws. It is a set of partial differential equations associated with the relationships of electromagnetic field (**E**, **H**) distributing into space and varying into time, the distribution of currents and charges (**J**, ρ) and material properties (μ, σ) [1,2].

In order to compute the electromagnetic field (i.e., solutions to Maxwell's equations), an analytical method or a numerical method can be used. For these devices with simple geometries, a correct analytical solution can be identified. However, in some general cases, especially for the complex structure of electrical devices, the numerical methods are used as a sole solution. The numeric methods applied in modelling of electromagnetic fields can be divided into two categories: finite methods like FEM (Finite Element Method), FVM (Finite Volume Method) and the numerical integration methods such as BEM (Boundary Elements method), MoM (Methods of Moments), PEEC method (Partial Element Equivalent Circuit) [3–13]. The choice of an appropriate method totally depends on the physical phenomena that need modelling such as the physical phenomena in high or low frequencies, with or without magnetic material, considering the effect of inductance or capacitance, external excitation source. However, there is no universal, optimal method for these problems, and the choice of the most suitable method depends on the nature of the electrical devices and their operation range.

Generally, the electromagnetic modelling of structures including thin shell model and thin wires is a complex problem in the fields of electrical engineering. Its geometry is characterized by a high ratio of the length and the thickness. Thus, the use of a volume mesh leads to a large number of elements and/that makes this method expensive and time consuming to apply to practical devices. Furthermore, due to the skin effect, when the skin depth δis thinner than the thickness *e*, the size of the studied meshes must be smaller than that of the skin depth. This leads to the difficulties when the thickness of some parts of the structures is very small in comparison with the overall size of the actual devices.

The shell element formulations have recently been developed by many authors in order to cope with difficulties in thin plate models, e.g., with the boundary element method (BEM) [3], with the finite element method (FEM) [5] and with the integral methods [11,14–16]. In [11], we presented a coupling method between the Partial Element Equivalent Circuit (PEEC) and an integro-differential method to model the devices that include thin electromagnetic shells and complicated conductor systems. However, this coupling formulation cannot be applied to problems when the skin depth is low in comparison with the thickness of the thin shell. In [14], the authors presented a hybrid of volume and surface integral formulations for the eddy current solution of the conductive regions with arbitrary geometry. However, this formulation cannot be used for the magnetic material and as in [11], the skin effect is not yet considered. In [16], we developed a general shell element formulation for modelling of thin magnetic and conductive regions. The thin shell is modelled by an integral method to avoid meshing the air region. As in [3,5,15], the field variation through the thickness is considered. Based on a simple discretization of the shell averaged surface, the number of unknowns is greatly reduced. However, this formulation is not suitable for model systems with complex conductive structures.

In this paper, a hybrid integral formulation is proposed to allow the modelling of an inhomogeneous structure constituted by conductors and thin magnetic and conductive shells in the general case (δ > *e* or δ ≈ *e* or δ < *e*). The modelling of thin magnetic and conductive shell regions is determined thanks to the integral formulation in [15,16]. A method that allows the modelling of the contributions of the inductors fed with the alternating currents will be presented in Section 3. The coupling formulations for integrating the interaction between conductors and thin shells material will be fully developed in Section 4. Finally, two numerical examples are presented in Section 5. The results obtained from this formulation are compared with those obtained from the FEM. Strong and weak points of our formulation are also analyzed.

## **2. Thin Shell Equation**

A thin electromagnetic shell (Ω) with thickness *e*, conductivity σ and linear permeability μ<sup>r</sup> in this study is illustrated in Figure 1, where Г<sup>1</sup> and Г<sup>2</sup> are the surface boundaries of the shell with the air region.

**Figure 1.** Thin region and associated notations [16].

The electromagnetic behavior for the side "1" of the shell regions is represented by [15,16]:

$$\int\_{\Gamma\_1} \mathbf{grad}\_s w \cdot (a\mathbf{H}\_{1s} - \beta \mathbf{H}\_{2s}) \mathrm{d}\Gamma + j\omega \int\_{\Gamma\_1} w \cdot \mathbf{B}\_1 \cdot \mathbf{n}\_1 \mathrm{d}\Gamma = 0,\tag{1}$$

where δ is the skin depth associated to the shell; *a* = (1 + *j*)/δ; α = *a*/σth(*ae*), β = *a*/σsh(*ae*); *w* is a set of nodal surface weighting functions; **H**1s and **H**2s are the tangential magnetic fields on both sides Г1, Г<sup>2</sup> and **n**<sup>1</sup> is the normal vector of the side Г1.

The symmetrical equation corresponding to the other side of the shell is obtained by simply exchanging the indices "1" and "2":

$$\int\_{\Gamma\_2} \mathbf{grad}\_{\mathbf{s}} w \cdot (\alpha \mathbf{H}\_{2\mathbf{s}} - \beta \mathbf{H}\_{1\mathbf{s}}) \mathrm{d}\Gamma + j\omega \int\_{\Gamma\_2} w \cdot \mathbf{B}\_2 \cdot \mathbf{n}\_2 \mathrm{d}\Gamma = 0,\tag{2}$$

where **n2** is the normal vector corresponding the side "2" of the shell.

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

Equations (1) and (2) are written on the averaged surface Г of the thin shell region. Subtracting Equations (2) and (1) leads to:

$$\delta(\boldsymbol{\alpha} + \boldsymbol{\beta}) \int\_{\Gamma} \mathbf{grad}\_{\boldsymbol{s}} w(\mathbf{H}\_{2s} - \mathbf{H}\_{1s}) \mathrm{d}\Gamma + j\boldsymbol{\omega} \int\_{\Gamma} w(\mathbf{2} \cdot \mathbf{B}\_{\mathbf{s}} \cdot \mathbf{n}) \mathrm{d}\Gamma = 0,\tag{3}$$

where **B**<sup>a</sup> = (**B**<sup>1</sup> + **B**2)/2 is defined as the averaged induction and **n** denotes the normal vector of the surface Г (Figure 1).

The expressions related to the tangential magnetic fields on both sides and the outside of the shell are:

$$\mathbf{H}\_{1s} = \mathbf{H}\_{01s} - \mathbf{grad}\_s \phi\_1; \ \mathbf{H}\_{2s} = \mathbf{H}\_{02s} - \mathbf{grad}\_s \phi\_2 \tag{4}$$

where **H**01s and **H**02s are the tangential fields generated by the inductors at side Г1, Г<sup>2</sup> and φ1, φ<sup>2</sup> denote the corresponding reduced scalar magnetic potentials.

By using (4) and assuming the small variations of **H**0s through the thickness of the thin shell, Equation (3) becomes:

$$(\alpha + \beta) \int\_{\Gamma} \mathbf{grad}\_{\mathsf{S}} w \cdot \mathbf{grad}\_{\mathsf{S}} \Delta \phi \mathrm{d}\Gamma + 2j\omega \int\_{\Gamma} w \cdot \mathbf{B}\_{\mathsf{a}} \cdot \mathbf{nd} \Gamma = 0,\tag{5}$$

where Δφ = φ<sup>1</sup> − φ<sup>2</sup> is the scalar magnetic potential discontinuity.

Using magnetization law of the linear material, equation (5) can be rewritten as:

$$\delta \left( a + \beta \right) \int\_{\Gamma} \mathbf{grad}\_{\delta} w \cdot \mathbf{grad}\_{\delta} \Delta \phi \mathrm{d}\Gamma + 2j\omega \frac{\mu\_{0}\mu\_{r}}{\mu\_{r} - 1} \int\_{\Gamma} w \cdot \mathbf{M}\_{\mathfrak{d}} \cdot \mathbf{nd} \Gamma = 0,\tag{6}$$

where **M**<sup>a</sup> = (**M**<sup>1</sup> + **M**2)/2 is defined as the averaged magnetization.

On the averaged surface Г, the total magnetic field **H**<sup>a</sup> is the sum of the inductor field **H**0, the field created by magnetization **H**<sup>M</sup> and the field created by the eddy currents flowing in the shell **H**EC:

$$\mathbf{H}\_{\mathbf{4}}(P) = \mathbf{H}\_{0}(P) + \mathbf{H}\_{\mathbf{M}}(P) + \mathbf{H}\_{\text{EC}}(P) \tag{7}$$

where *P* is a point located in the surface region Г.

The field **H**<sup>M</sup> is equal to [10]:

$$\mathbf{H}\_{\mathbf{M}}(P) = -\mathbf{grad}\frac{1}{4\pi} \int\_{\Omega} \frac{(\mathbf{M} \cdot \mathbf{r})}{\mathbf{r}^3} d\Omega \tag{8}$$

where **r** is the vector between the integration point on Ω and the point *P*.

By using Biot and Savart law, the field **H**EC can be written as:

$$\mathbf{H}\_{\rm EC}(P) = \frac{1}{4\pi} \int\_{\Omega} \frac{\mathbf{J} \times \mathbf{r}}{\mathbf{r}^3} d\Omega,\tag{9}$$

where **J** denotes the eddy current density of the thin shell.

The integration of the tangent magnetization through the depth of the thin shell is written as follows [3,5,16]:

$$\int\_{-\epsilon/2}^{+\epsilon/2} \mathbf{M} \, \mathbf{dz} = \frac{\overline{\mathbf{G}}(\mathbf{M}\_1 + \mathbf{M}\_2)}{2} = \overline{\mathbf{G}} \mathbf{M}\_{\mathbf{3}} \tag{10}$$

where G = th[(1 + *j*)*e*/2δ]/[(1 + *j*)*e*/2δ].

Using (10), Equation (8) becomes:

$$\mathbf{H}\_{\rm M}(P) = -\mathbf{grad}\frac{\overline{\mathbf{G}}}{4\pi} \int\_{\Gamma} \frac{(\mathbf{M}\_{\rm a} \cdot \mathbf{r})}{\mathbf{r}^3} d\Gamma. \tag{11}$$

The equivalent shell current **K** is defined as [3]:

$$\mathbf{K} = \int\_{-\epsilon/2}^{+\epsilon/2} \mathbf{J} d\mathbf{z} = \mathbf{n} \times \mathbf{grad} \Delta \phi. \tag{12}$$

Using (12), Equation (9) becomes:

$$\mathbf{H}\_{\rm EC}(\mathbf{P}) = \frac{1}{4\pi} \int\_{\Gamma} \frac{\mathbf{n} \times \mathbf{grad}\Delta\phi \times \mathbf{r}}{\mathbf{r}^3} d\Gamma. \tag{13}$$

Combining (11) with (12), Equation (7) is rewritten as:

$$\mathbf{H}\_{\mathsf{a}}(P) = \mathbf{H}\_{0}(P) - \mathbf{grad}\frac{\overline{\mathbf{G}}}{4\pi} \int\_{\Gamma} \frac{(\mathbf{M}\_{\mathsf{a}} \cdot \mathbf{r})}{\mathbf{r}^{3}} d\Gamma + \frac{1}{4\pi} \int\_{\Gamma} \frac{\mathbf{n} \times \mathbf{grad}\Delta\phi \times \mathbf{r}}{\mathbf{r}^{3}} d\Gamma. \tag{14}$$

Let us consider a linear magnetic law for the material **M**a(*P*) = (μ*<sup>r</sup>* − 1)**H**a(*P*), then Equation (14) becomes [10,11]:

$$\frac{\mathbf{M}\_{\mathbf{a}}(P)}{\mu\_{r} - 1} = \mathbf{H}\_{0}(\mathbf{P}) - \mathbf{grad}\frac{\overline{\mathbf{G}}}{4\pi} \int\_{\Gamma} \frac{(\mathbf{M}\_{\mathbf{a}} \cdot \mathbf{r})}{\mathbf{r}^{3}} d\Gamma + \frac{1}{4\pi} \int\_{\Gamma} \frac{\mathbf{n} \times \mathbf{grad}\Delta\phi \times \mathbf{r}}{\mathbf{r}^{3}} d\Gamma,\tag{15}$$

In [16], Equations (6) and (15) have been resolved and validated by the comparison with the axisymmetric FEM and the FEM 3D with the shell elements. Let us note that in the proposed equations, the integrals are determined only by the surface numerical integration. Therefore, the 3D implementation is simple and reliable. The comparison results with different ratios (*e*/δ) have demonstrated the ability and the advantages of the method.

#### **3. Conductor System Modelling**

We now consider *m* volume conductors fed with alternating sources placed in an air region. The external electric field incident at point *P* can be written [7–9]:

$$\mathbf{E}\_{\rm ext}(P) = \frac{\mathbf{J}\_c(P)}{\sigma} + j\omega \mathbf{A}(P) + \mathbf{grad}\mathbf{V}(P),\tag{16}$$

where **A**(*P*) is the vector potential and V(*P*) is the scalar potential.

The magnetic vector potential generated by the current density **Jc** is:

$$\mathbf{A}(P) = \sum\_{k=1}^{m} \frac{\mu\_0}{4\pi} \int\_{\Omega c\_k} \frac{\mathbf{J}\_c}{\mathbf{r}} d\Omega,\tag{17}$$

where Ω*ck* is the volume of the conductor *k*; r is the distance between the integration point on Ω*ck* and the point *P*.

In the context of the quasi-stationary regime operated at the frequency range up to ten MHz, it is possible to neglect the capacitive effect. Equation (16) is written in the form:

$$\mathbf{E}\_{\rm ext}(P) = \frac{\mathbf{J}\_{\rm c}(P)}{\sigma} + j\omega \sum\_{k=1}^{m} \frac{\mu\_{0}}{4\pi} \int\_{\Omega \mathbf{x}\_{k}} \frac{\mathbf{J}\_{\rm c}}{\mathbf{r}} d\Omega. \tag{18}$$

The PEEC method is particularly pertinent and reliable to solve this kind of problem. Let us assume that the current density in each conductor is uniform. For each conductor, Equation (18) can be associated to the electrical equivalent circuit presented by the self and mutual inductances in series with the resistances [8,9]. In most respects, what we know about this approach is difficult to the magnetic material or the magnetic conductive material. However, the mesh of the air region can be neglected and some conventional SPICEs-like or Saber circuit solvers can be employed to analyze the equivalent circuit. As a result, the behavior of several electrical devices that have electrical interconnections can be studied in only one system. With all of the advantages described above, the PEEC is well-suited for modelling real industrial devices [17–21].

Despite these clear advantages, this technique has a major drawback due to the necessity to store a fully dense matrix inherent to the use of the integral formulation in Maxwell's equation. Consequently, this approach is strongly limited in modelling large-scale electrical devices requiring substantial meshes. Let us denote N as the number of unknown then the needed memory for the fully dense matrix storage leads to O(N2) and the computational costs of a direct solver (LU-decomposition) increases in proportion to O(N3). For example, a large size modelling problem containing 50,000 unknowns requires at least 50,000 × 50,000 × 8 bytes or approximately 20 GB in the memory to store a fully dense matrix and the computational costs of finding one solution by direct solver are roughly several weeks. In order to study the characteristic of the devices in a frequency range or in a period of time by the PEEC method, the total computational costs must be multiplied with the number of frequency or the time steps. Moreover, none of the circuit solvers can handle the equivalent circuit composing of 50,000 × 50,000 basic circuit elements RLMC (resistor, inductor, partial inductor, and capacitor). Overall, in both cases the computational time tends to infinity, or it is impossible to solve. To solve this problem, the algorithms coupling different matrix compression algorithms or the model order reduction technique with integral methods have been developed with the purpose of limiting matrix storage [19,21].

In addition, the expansion of the PEEC method for more general problems has been investigated [11,18,20]. However, this method is still difficult for considering special structures such as a circular coil, thick circular coil, and thin disk coil. This drawback can be easily improved by using the semi-analytic methods developed in [22,23].
