3-D Integral Formulation for Thin Electromagnetic Shells Coupled with an External Circuit

The aim of this article is to present a hybrid integral formulation for modelling structures made by conductors and thin electromagnetic shell models. Based on the principle of shell elements, the proposed method provides a solution to various problems without meshing the air regions, and at the same time helps to take care of the skin effect. By integrating the system of circuit equations, the method presented in this paper can also model the conductor structures. In addition, the equations describing the interaction between the conductors and the thin shell are also developed. Finally, the formulation is validated via an axisymmetric finite element method and the obtained results are compared with those implemented from another shell formulation.


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][4][5][6][7][8][9][10][11][12][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][15][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.

Thin Shell Equation
A thin electromagnetic shell (Ω) with thickness e, conductivity σ and linear permeability µ r in this study is illustrated in Figure 1, where Г 1 and Г 2 are the surface boundaries of the shell with the air region.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 2 of 13 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][15][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.

Thin Shell Equation
A thin electromagnetic shell (Ω) with thickness e, conductivity σ and linear permeability μr in this study is illustrated in Figure 1, where Г1 and Г2 are the surface boundaries of the shell with the air region. The electromagnetic behavior for the side "1" of the shell regions is represented by [15,16]: where is the skin depth associated to the shell; = (1 + )/ ; = / th( ), = / sh( ); w is a set of nodal surface weighting functions; H1s and H2s are the tangential magnetic fields on both sides Г1, Г2 and n1 is the normal vector of the side Г1. The electromagnetic behavior for the side "1" of the shell regions is represented by [15,16]: 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 , Г 2 and n 1 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": where n 2 is the normal vector corresponding the side "2" of the shell.
Equations (1) and (2) are written on the averaged surface Г of the thin shell region. Subtracting Equations (2) and (1) leads to: where B a = (B 1 + 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: where H 01s and H 02s are the tangential fields generated by the inductors at side Г 1 , Г 2 and φ 1 , φ 2 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: where ∆φ = φ 1 − φ 2 is the scalar magnetic potential discontinuity. Using magnetization law of the linear material, equation (5) can be rewritten as: where M a = (M 1 + M 2 )/2 is defined as the averaged magnetization. On the averaged surface Г, the total magnetic field H a is the sum of the inductor field H 0 , the field created by magnetization H M and the field created by the eddy currents flowing in the shell H EC : where P is a point located in the surface region Г.
The field H M is equal to [10]: 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: 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]: where Using (10), Equation (8) becomes: The equivalent shell current K is defined as [3]: Using (12), Equation (9) becomes: Combining (11) with (12), Equation (7) is rewritten as: Let us consider a linear magnetic law for the material M a (P) = (µ r − 1)H a (P), then Equation (14) becomes [10,11]: 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.

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][8][9]: where A(P) is the vector potential and V(P) is the scalar potential. The magnetic vector potential generated by the current density J c is: Appl. Sci. 2020, 10, 4284 where Ωc k is the volume of the conductor k; r is the distance between the integration point on Ωc k 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: 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][18][19][20][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].

Coupling Thin Shell with Circuit Equation
We now consider a general system composed of m conductors and thin magnetic conductive shell (volume Ω and average surface Γ) (Figure 2). Let us assume that the current density in each conductor is uniform.

Coupling Thin Shell with Circuit Equation
We now consider a general system composed of m conductors and thin magnetic conductive shell (volume Ω and average surface Γ) (Figure 2). Let us assume that the current density in each conductor is uniform.

Influence of the Conductor Current on the Thin Region
The term H0 in Equation (15) can be computed by using Biot and Savart law: where Ik and Sk are respectively the current and the cross section of the conductor k; Ωc is the volume of conductor k; is the vector unit of current direction in the conductor k and r is the vector between the integration point on Ωc and the point P where the field is expressed.
Equation (15) is then rewritten as:

Influence of the Conductor Current on the Thin Region
The term H 0 in Equation (15) can be computed by using Biot and Savart law: where I k and S k are respectively the current and the cross section of the conductor k; Ωc k is the volume of conductor k; l k is the vector unit of current direction in the conductor k and r is the vector between the integration point on Ωc k and the point P where the field is expressed. Equation (15) is then rewritten as:

Influence of Thin Shell Magnetization on the Conductor
In order to take into account the influence of the field created by the thin shell magnetization, we have to integrate the magnetic vector potential A M generated by the magnetized thin shell on the conductor. This magnetic vector potential is expressed as [11,24]: where r denotes the vector between the integration point and the point P.

Influence of Thin Shell Eddy Current on the Conductor
The eddy current J in the thin shell generates a magnetic potential vector on the conductor: where r is the distance between the integration point on Ω and the point P.

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.
The Galerkine projection method is then applied to Equation (6), and we get the following matrix system [11,15,16]: where: [M] is a complex vector of dimension 3n; [∆Φ] is a complex vector of dimension p; [A] is a (p × 3n) matrix expressed as follows: [B] is a (p × p) sparse matrix and can be written as: 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]: where [,] is defined as the scalar product operator and u i is the vector basis of the element i and r i is the vector between the integration point on Γ j to the centroid of the element i.
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: where R k is the resistance of the k th conductor; m ik is the mutual inductance between the two conductors k, i; m f ik is defined as the mutual inductance between the magnetization M ai of the shell element and the conductor k; m q 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: where u i is a basis vector of magnetizations; n i is the normal vector of the element Γ i and l k is the vector unit of the current direction in the conductor k. Let us note that the mutual inductance terms m ik can be calculated by the well-known semi-analytic formulas in [7][8][9]22,23]. Whereas, the terms m f ik , m q 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: where Finally, the algebraic linear systems (26), (29) and (37) is considered as p + 3n + m complex unknowns: µ r −1 + [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 ]/(µ r − 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.

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.

Validation Through an Academic Example
In the first example, a thin magnetic conductive disk is considered with the circular conductor fed by a voltage of 1V and a frequency of 10Hz (Figure 3).
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.

Validation Through an Academic Example
In the first example, a thin magnetic conductive disk is considered with the circular conductor fed by a voltage of 1V and a frequency of 10Hz (Figure 3). This academic example is solved by three numerical methods, where the first one is the axisymmetric FEM, the second one is a shell element formulation implemented in 3D FEM method [5] and the last one is the proposed integral method with the surface elements. In order to valid our method and to compare different approaches, we mainly focus on the computed current in the conductor and the magnetic field in the air region close to the device (calculated on the path AB, for A (0.25; 0; 0.1) and B (0.25; 0; 0.25)).
Let us note that the problem must be meshed very finely to have an accurate result with the FEM 3D (Table 1 and Figure 4) because of the high variations of the fields around the conductor. The current values greatly vary according to the number of the elements.
If the axisymmetric FEM is considered as our reference, the coupling integral method leads to an error of 0.1% (Table 1). Our method leads to more accurate results than the same shell element formulation but considering the air region treated with the 3D FEM.  [5] and the last one is the proposed integral method with the surface elements. In order to valid our method and to compare different approaches, we mainly focus on the computed current in the conductor and the magnetic field in the air region close to the device (calculated on the path AB, for A (0.25; 0; 0.1) and B (0.25; 0; 0.25)). Let us note that the problem must be meshed very finely to have an accurate result with the FEM 3D (Table 1 and Figure 4) because of the high variations of the fields around the conductor. The current values greatly vary according to the number of the elements.  If the axisymmetric FEM is considered as our reference, the coupling integral method leads to an error of 0.1% (Table 1). Our method leads to more accurate results than the same shell element formulation but considering the air region treated with the 3D FEM.

A Pratical Device Example
The second test problem is the modelling of a practical device proposed by the EDF (Electricité de France) [26]. The power station "Folies" is equipped with a three-phase reactance to limit shortcircuit currents. In this case, the current in each reactance is 1000 A phase-shifted with 120 degrees, the frequency is 50 Hz ( Figure 5).

A Pratical Device Example
The second test problem is the modelling of a practical device proposed by the EDF (Electricité de France) [26]. The power station "Folies" is equipped with a three-phase reactance to limit short-circuit currents. In this case, the current in each reactance is 1000 A phase-shifted with 120 degrees, the frequency is 50 Hz ( Figure 5).  In nominal conditions, a three-phase reactance generate a leakage induction in the neighborhoods. In order to minimize this electromagnetic disturbance, an electromagnetic shielding and a passive loop are added between the magnetic field source and the protected area (Figure 7). The device parameters are presented as below: In nominal conditions, a three-phase reactance generate a leakage induction in the neighborhoods. In order to minimize this electromagnetic disturbance, an electromagnetic shielding and a passive loop are added between the magnetic field source and the protected area (Figure 7). The device parameters are presented as below: -Passive loop: Permeability µ r = 1, Electrical conductivity σ = 3.03 × 10 7 S/m Section radius r s = 9.25 mm Z loop = 3.85 m -Shielding: Permeability µr = 20,000, Electrical conductivity σ = 2.2 × 106 S/m Thickness = 3.5 mm Zshielding = 5.15 m The skin depth δ = 0.339mm is thinner than the shielding's thickness.
The last case is tested by our integral method and the FEM 3D with shell elements. Let us note that this example is in 3-dimensional space and cannot be modelled by the 2D FEM. The current values in the passive loop and the current distribution in the shells are also compared.
The obtained results from the coupling method converge quite close to the current values presented in Table 2. Figure 8 also shows that the surface distribution of the current in the shell is quite similar to the two methods. The results achieved by the coupling integral method are very encouraging. The convergence is reached with few elements (about 1000 elements). The skin depth δ = 0.339mm is thinner than the shielding's thickness.  The last case is tested by our integral method and the FEM 3D with shell elements. Let us note that this example is in 3-dimensional space and cannot be modelled by the 2D FEM. The current values in the passive loop and the current distribution in the shells are also compared.
The obtained results from the coupling method converge quite close to the current values presented in Table 2. Figure 8 also shows that the surface distribution of the current in the shell is quite similar to the two methods. The results achieved by the coupling integral method are very encouraging. The convergence is reached with few elements (about 1000 elements).   o The skin depth δ = 0.339mm is thinner than the shielding's thickness.  The last case is tested by our integral method and the FEM 3D with shell elements. Let us note that this example is in 3-dimensional space and cannot be modelled by the 2D FEM. The current values in the passive loop and the current distribution in the shells are also compared.
The obtained results from the coupling method converge quite close to the current values presented in Table 2. Figure 8 also shows that the surface distribution of the current in the shell is quite similar to the two methods. The results achieved by the coupling integral method are very encouraging. The convergence is reached with few elements (about 1000 elements).   However, some matrices of Equation (38) are fully populated and compression algorithms must be applied if there is a large number of elements. The coupling of the model order reduction techniques or the matrix compression algorithms with some integral methods clearly demonstrates its efficiency. For example, the coupling of a matrix compression algorithm like the FMM and the MoM or the PEEC method has reduced the computation time and the memory requirements down to more than 10 times and the compressed ratio is more than 80 percent [18,27]. Moreover, the acquired model can be reused to build a real circuit, which is easy to employ in all conventional SPICEs-like circuit solvers [28,29].

Conclusions
In this paper, we have presented a coupling integral method in order to model thin magnetic and conductive regions which can be coupled with an external electric circuit. This coupling enables the model of conductors with complex shapes, and various skin effects across the thickness of the thin shell are taken into account. Two numerical examples have been presented and the results highlighted the accuracy of the solution provided by our proposed method. In our suggestions, the problem dense matrix and the full memory need could be solved with the compression algorithms.

Conflicts of Interest:
The authors declare that there are no conflicts of interest regarding the publication of this paper.

Abbreviations
The following abbreviations are used in this paper: