A FEM-Green Approach for Magnetic Field Problems with Open Boundaries

: A new ﬁnite element method/boundary element method (FEM/BEM) scheme is proposed for the solution of the 2D magnetic static and quasi-static problems with unbounded domains. The novelty is an original approach in the treatment of the outer region. The related domain integral is eliminated at the discrete level by using the ﬁnite element approximation of the fundamental solutions (Green’s functions) at every node of the related mesh. This “FEM-Green” approach replaces the standard boundary element method. It is simpler to implement because no integration on the boundary of the domain is required. Then, the method leads to a substantially reduced computational burden. Moreover, the coupling with ﬁnite elements is more natural since it is based on the same Galerkin approximation. Some examples with open boundary and nonlinear materials are presented and compared with the standard ﬁnite element method.


Introduction
A large number of electromagnetic problems involve finding fields in the vicinity of devices embedded in an infinite domain. However, classical volume discretization methods, such as the finite element method, are inherently inapplicable. Various strategies have been developed over the years, such as simple truncation of outer boundaries, ballooning, conformal mappings, infinite elements, hybrid finite element method (FEM)-boundary element method (BEM) schemes, and others. An excellent review is given in [1]. More specifically, FEM/BEM methods exploit the flexibility of the FEM for the treatment of nonhomogeneous and nonlinear materials, and the ability of the BEM to consider open regions [2]. The FEM is generally based on the Galerkin approximation of a weak or variational formulation, while the standard BEM is normally deduced from the Green's second identity for reducing domain integrals into boundary ones [3].
Instead of this classical BEM approach, in this paper, we propose an alternative socalled FEM-Green formulation as first described by the author in [4,5]. It is obtained by using the FEM approximation of fundamental solutions (Green's functions) [6,7] and consists of a discrete equivalent of the second Green's identity. More precisely, the domain integral of the open domain is naturally canceled by a combination of the Galerkin formulation of the boundary value problem under study, and the one associated with the fundamental solution for the Dirac delta loading every node of a finite element mesh of this region. The layer of finite elements that are adjacent to the boundary of the FEM region is only required with our method. Any extended mesh in some parts of the outer domain is used for the computation of fields at a postprocessing step.
The implementation of this new scheme is simpler as no boundary integration, either analytical or numerical (Gaussian quadrature), on the FEM/BEM interface is required and the coupling with the standard FEM used in the inner part of the problem appears to be extremely natural. In the two previous papers of the author, the mathematical developments 2 of 11 were neither presented with a source term nor in the context of magnetic field problems. Indeed, those problems have their own specificities such as a large difference of permeability between air and ferromagnetic materials. Moreover, time-harmonic eddy-current problems are described here by a complex-valued equation, which was not considered before. Thus, it is worth exploring the applicability of our method to such problems. It is the objective of the present paper. Therefore, the method is applied to 2D magnetostatic and eddy-current open boundary problems and some practical examples with open boundary and nonlinear domains are presented in the last section in order to assess the performance of our method.

Mathematical Model Description
Various FEM/BEM hybrid methods are used for the computation of open boundary problems where the classical BEM formulation is used for the treatment of the open region. Green's second identity is normally employed to derive a boundary-only formulation, avoiding the discretization of the domain. The principle of the proposed method is different and is based upon a special FEM treatment using the fundamental solution (Green's function) of the boundary value problem. This allows for a direct elimination of the domain integral at the discrete level through a Galerkin approach.
The method is presented in the context of the modified Helmholtz problem describing time-harmonic two-dimensional quasi-static problems by using magnetic vector potential (MVP) A and the magnetostatics situation as a particular case. Magnetic vector potential possesses only a longitudinally directed component denoted as A (also commonly called magnetic vector potential), so that the magnetostatic vector Poisson equation degenerates to its scalar counterpart that is stated here below. Since the 2D aspect is not mandatory for the method, it could be extended to 3D problems. This speculation is valid for scalar potential problems, e.g., in electrostatics using the scalar electric potential, as it is mostly the case, or in 3D magnetostatics using a magnetic scalar potential [8]. However, it should be fully analyzed for magnetic field problems using full MVP. Figure 1 depicts the general 2D configuration where Ω 1 is the open region (air), including the source region Ω s with current density J s , and Ω 2 is the inner region that can be nonhomogeneous and/or nonlinear and/or with eddy-currents, depending on whether a magnetostatic or an eddy-current analysis is considered.
Permeability µ and conductivity σ are assumed as essential material properties and we consider no Dirichlet or Neumann boundary condition for the sake of simplicity in the mathematical description of the method. The governing equations are: and while the following continuity conditions hold on the interface Γ int : Potential A and current density J s are represented in phasor form with angular frequency ω in case of a time-harmonic analysis.

The FEM-Green Formulation
The FEM-Green formulation is applied to the open region Ω 1 and is derived from a finite element mesh. It is first assumed that domain Ω 1 is bounded by the interface Γ int and a remote boundary Γ ∞ where the homogeneous Dirichlet boundary condition (A = 0) is considered (Figure 1). The mathematical treatment will show that the domain integral vanishes so that only a layer of finite elements along the interface Γ int is required and there will be no problem when outer boundaryΓ ∞ extends to infinity.
for scalar potential problems, e.g., in electrostatics using the scalar electric potential, as it is mostly the case, or in 3D magnetostatics using a magnetic scalar potential [8]. However, it should be fully analyzed for magnetic field problems using full MVP. Figure 1 depicts the general 2D configuration where Ω is the open region (air), including the source region Ω with current density , and Ω 2 is the inner region that can be nonhomogeneous and/or nonlinear and/or with eddy-currents, depending on whether a magnetostatic or an eddy-current analysis is considered. Mathematical treatment used here is classical and can be found in many standard books such as [9]. Consider the well-known Galerkin problem related to Poisson's Equation (1a): where the N j 's are the classical interpolation functions defined on finite element mesheŝ Ω 1 andΩ s of the domains Ω 1 and Ω s , respectively ( Figure 1). We will assume firstorder triangular elements for the sake of clarity in the presentation of the method. The approximate solutionÂ is given by: The combination of Equations (3) and (4) leads to a sparse system of linear equations, of which the matrix entries are of the form: Let us denote It is demonstrated in [10] that this term corresponds physically to the inward flux of ∇Â across the "box" Σ j associated with node j in the dual mesh obtained from the barycentric subdivision of the primal 2D mesh. In particular, at a node j belonging to the discretization of the interface Γ int between Ω 1 and Ω 2 , the quantity ΦÂ ,Σ j is interpreted as the flux through the "cap" Σ j , as shown in Figure 2.
However, this interpretation is only valid for first-order elements, and it is not the case of higher-order elements where it is a mere mathematical equivalent variable. Yet, this variable obeys a continuity condition across the interface, as will be shown in Section 2.2. Parameters ΦÂ ,Σ j are used here instead of the approximation of the normal derivative that is normally present in the classical boundary element method [3].
Mathematics 2021, 9, x FOR PEER REVIEW 4 of 12 centric subdivision of the primal 2D mesh. In particular, at a node belonging to the discretization of the interface Γ between Ω and Ω , the quantity Φ , is interpreted as the flux through the "cap" Σ , as shown in Figure 2. However, this interpretation is only valid for first-order elements, and it is not the case of higher-order elements where it is a mere mathematical equivalent variable. Yet, this variable obeys a continuity condition across the interface, as will be shown in Section 2.2. Parameters Φ , are used here instead of the approximation of the normal derivative that is normally present in the classical boundary element method [3].
As in the standard BEM, the 2D fundamental solution = −1/2 ln of Laplace's equation is exploited in order to eliminate the domain integral. Let us recall that it is obtained from the equation where is the Dirac delta function at any point ( Figure 1). A FEM solution exists [6,7], and is readily derived from the Galerkin problem associated with the formulation (7): where the FEM solution is interpolated as An illustrative example of the FEM solution of a Green's function on an arbitrary domain Ω is given in Figure 3.  As in the standard BEM, the 2D fundamental solution G i = −1/2π ln r i of Laplace's equation is exploited in order to eliminate the domain integral. Let us recall that it is obtained from the equation where δ i is the Dirac delta function at any point i ( Figure 1). A FEM solution exists [6,7], and is readily derived from the Galerkin problem associated with the formulation (7): where the FEM solutionĜ i is interpolated aŝ An illustrative example of the FEM solutionĜ i of a Green's function G i on an arbitrary domain Ω is given in Figure 3. However, this interpretation is only valid for first-order elements, and it is not the case of higher-order elements where it is a mere mathematical equivalent variable. Yet, this variable obeys a continuity condition across the interface, as will be shown in Section 2.2. Parameters Φ , are used here instead of the approximation of the normal derivative that is normally present in the classical boundary element method [3].
As in the standard BEM, the 2D fundamental solution = −1/2 ln of Laplace's equation is exploited in order to eliminate the domain integral. Let us recall that it is obtained from the equation where is the Dirac delta function at any point ( Figure 1). A FEM solution exists [6,7], and is readily derived from the Galerkin problem associated with the formulation (7): where the FEM solution is interpolated as An illustrative example of the FEM solution of a Green's function on an arbitrary domain Ω is given in Figure 3.  In the same way as forÂ given by (6), we can write the nodal approximation: Note that ΦĜ i ,Σ i = 1, as a consequence of Gauss's law for any internal node i of the meshΩ 1 . If we combine Equations (3)-(6) and (9), a little algebra suffices to show that: Mathematics 2021, 9, 1662 5 of 11 By performing integration of the source term, (11) becomes: where T s,j is the set of triangles ∆ j of the source regionΩ s sharing the node j and ∆ j denotes the area of the triangle. Conversely, using (4), adapted toĜ i and Equation (10) instead of (6), we can easily write: Now, by equating (12) to (13), we can write the following expression: where the domain contribution is eliminated as announced previously, and∂ Ω 1 has been replaced byΓ int since it is now licit to extendΓ ∞ the outer boundary at infinity. Equation (14) can be regarded as a discrete equivalent of Green's second identity that has been derived directly from the Galerkin approach. Henceforth, it is necessary to write (14) for all the nodes i belonging to the boundaryΓ int , in order to derive a consistent linear system of equations. However, due to the singularity of the G i function, the boundary condition of the problem (7) must be changed as: so that the Galerkin formulation (8) is replaced by: where δ ij is the common Kronecker symbol and c i is the classical geometric factor of standard BEM, i.e., it is equal to the ratio θ i /2π, where θ i is the internal (with respect to the FEM region Ω 2 ) angle at node i [3]. It is easy to show that, in the FEM-Green context: that can again be interpreted as Gauss's law at the discrete level. Finally, the FEM-Green scheme leads to the system of the n simultaneous equations: where potentialÂ j and flux ΦÂ ,Σ j values are the unknown parameters along the interfaceΓ int .
At this point, the solution of (18) requires the computation of the FEM approximationŝ G i that is highly time consuming and does not make sense sinceΩ 1 is unbounded in the context of open boundary problems. Then, coefficients ΦĜ i ,Σ j andĜ i are modified by using the exact fundamental solutions G i instead ofĜ i . More precisely the following interpolant should be considered: where r ij is the distance between nodes i and j. However, the infinite value G i (i) must be replaced by a value G * i,i that can be derived from the discrete Gauss's law (17), i.e., and:ĉ Parameterĉ i can be considered as an estimate of the geometric factor c i described above. By eliminating this parameter between Equations (20a) and (20b), after some algebra that we skip for the sake of conciseness, we obtain the expression of G * i,i : As outlined above, meshing the whole domainΩ 1 is not a necessity since the internal nodes are not involved in (18). In fact, that would not make sense since this outer region is unbounded. A single layer of finite elements along the boundary ∂Ω 1 (in gray in Figure 1) is sufficient. Any internal mesh of Ω 1 is used for field calculation at a postprocessing step.
By comparing with the implementation of the boundary element method, no cumbersome analytical or numerical (Gaussian quadrature) integration is required to compute the coefficients Φ G i,I ,Σ j , G i (j) and G * i,i of (18) so that an obvioussignificant reduction of the computational burden is expected. However, the computational effort to build the linear system (18) still scales as O n 2 as in classical BEM.
Lastly, note that the method can be applied to axisymmetrical problems where the fundamental solution is based on a complete elliptic integral of the first kind as it was shown by the author in [5].

FEM/FEM-Green Coupling
A complete set of equations associated with the hybrid FEM/FEM-Green is necessary to solve the whole problem. Then, the Galerkin problem related to the finite element domain Ω 2 has now to be derived. By referring again to Figure 1 and Equation (1b), the governing equation is either: for general nonlinear magnetostatic problems, or in case of time-harmonic eddy-current problems. The FEM problem is given by the respective Galerkin formulations The Galerkin equation written for a node j belonging to the interfaceΓ int in a conventional FEM method would be: The right-hand side of (25) vanishes for nodes j that are not adjacent to the source regionΩ s . By introducing the expression (6) of ΦÂ ,Σ j , Equation (25) becomes: Finally, a global system of algebraic equations of the whole problem is obtained by an assembling procedure of Equations (18) and (24), or (25), and (27). It may be expressed as the partitioned matrix form: where the unknown vector Φ i refers to the nodal flux values onΓ int , and vectorsû i andû 2 are related to the nodal potential values onΓ int andΩ 2 , respectively. Submatrices H and G represent the FEM-Green equations with the entries Φ G i,I ,Σ j and G i (j), or G * i,i , respectively. The submatrix S comes from the FEM contribution and 1 i is a unit matrix induced by (27). As it is the case for most FEM/BEM coupling methods, the global matrix has no particular structure, i.e., it is neither symmetric, nor positive definite. However, the G submatrices are symmetric. A general solver must be used for the solution of the system (28), but the optimization of this specific point has not been investigated in the paper. Finally, vector b relates to the source excitation J s that appears in the right-hand side of (18) and (27).

Numerical Results
The method has been applied to two magnetic problems to estimate the performance of our numerical scheme. In the first example, a comparison with the classical FEM/BEM technique and the FEM with a large truncated outer boundary has been carried out to assess the accuracy of the method. As in previous papers of the author [4,5], all the algorithms were implemented in the MATLAB ® environment on a standard desktop computer. The LU-decomposition was used for the solution of the involved linear systems. COMSOL Multiphysics ® software has been used for the generation of the various meshes required by the simulations and also in order to confirm the results obtained by the FEM.

Example 1-A Nonlinear Magnetostatic Problem
The first example is a 2D nonlinear magnetostatic problem defined on a C-shaped magnetic circuit Ω 2 with an excitation coil Ω s (J s = 10 7 A/m 2 ) in air region Ω 1 , as depicted in Figure 4.
As an illustration, a magnetic field plot obtained by using COMSOL Multiphysics ® is given in Figure 5. As classical FEM is used, the problem is encased in a circular box with a radius equal to about 15 times the mean size of the device. Magnetic induction along the line AB crossing the air gap, as shown in Figure 5, is plotted for the three methods in Figure 6. It is clear that FEM/BEM and FEM/FEM-Green methods are both in excellent agreement. Table 1 presents the magnetic energy restricted in an arbitrary rectangular domain of size 0.6 × 0.4 m centered on the origin of the x − y reference frame of the geometry and computed at a postprocessing step for both FEM/BEM and FEM/FEM-Green techniques. The number m of elements of magnetic and coil regions Ω 2 ∪ Ω s , is taken as a parameter. Convergence is observed and an extremely good agreement is obtained between the methods. A value of 74.53 J/m has been obtained with COMSOL in the same domain, by using a supplementary large external bounding rectangular box of size 10 × 10 m. This result is in coherence with our results. rithms were implemented in the MATLAB ® environment on a standard desktop computer. The LU-decomposition was used for the solution of the involved linear systems. COMSOL Multiphysics ® software has been used for the generation of the various meshes required by the simulations and also in order to confirm the results obtained by the FEM.

Example 1-A Nonlinear Magnetostatic Problem
The first example is a 2D nonlinear magnetostatic problem defined on a C-shaped magnetic circuit Ω with an excitation coil Ω ( = 10 A/m 2 ) in air region Ω , as depicted in Figure 4. As an illustration, a magnetic field plot obtained by using COMSOL Multiphysics ® is given in Figure 5. As classical FEM is used, the problem is encased in a circular box with a radius equal to about 15 times the mean size of the device. Magnetic induction along the line AB crossing the air gap, as shown in Figure 5, is plotted for the three methods in Figure 6. It is clear that FEM/BEM and FEM/FEM-Green methods are both in excellent agreement. Table 1 presents the magnetic energy restricted in an arbitrary rectangular domain of size 0.6 0.4 m centered on the origin of the − reference frame of the geometry and computed at a postprocessing step for both FEM/BEM and FEM/FEM-Green techniques. The number of elements of magnetic and coil regions Ω ∪ Ω , is taken as a parameter. Convergence is observed and an extremely good agreement is obtained between the methods. A value of 74.53 J/m has been obtained with COMSOL in the same domain, by using a supplementary large external bounding rectangular box of size 10 10 m. This result is in coherence with our results. As already emphasized in [5], there is a lower computational burden in the FEM-Green part of the system (28), compared with the equivalent BEM part. However, global CPU times measured for both methods for the proposed example are not vastly different due to the important part involved in the solutions of the system at every iteration.   As already emphasized in [5], there is a lower computational burden in the FEM-Green part of the system (28), compared with the equivalent BEM part. However, global CPU times measured for both methods for the proposed example are not vastly different due to the important part involved in the solutions of the system at every iteration.

Example 2-Eddy-Currents in a Conducting Magnetic Region
The second example deals with a conducting region Ω 2 made of copper (µ = 50 × µ 0 H/m and σ = 5.998 × 10 7 S/m) and driven by source currents (J s = 10 7 A/m 2 , frequency = 5 Hz) located in region Ω s , embedded in air region Ω 1 , as presented in Figure 7.

Example 2-Eddy-Currents in a Conducting Magnetic Region
The second example deals with a conducting region Ω made of copper ( = 50 H/m and = 5.998 10 S/m) and driven by source currents ( = 10 A/m 2 , frequency = 5 Hz) located in region Ω , embedded in air region Ω , as presented in Figure 7. Time-harmonic modified Helmholtz equation is applied to this problem. COMSOL Multiphysics ® solution is also given in Figure 7 where the induced current density is given by the statement: Real and imaginary parts of along the line AB crossing the conducting region is plotted in Figure 8 for FEM and FEM/FEM-Green methods (FEM-BEM is not applied here). Time-harmonic modified Helmholtz equation is applied to this problem. COMSOL Multiphysics ® solution is also given in Figure 7 where the induced current density is given by the statement: Real and imaginary parts of J i along the line AB crossing the conducting region is plotted in Figure 8 for FEM and FEM/FEM-Green methods (FEM-BEM is not applied here).
There is an excellent agreement between the results as the symbols • and * match the curves extremely well. The simulation has been realized with the number of triangular elements in conducting and source regions equal to 5811 as depicted in Figure 7. by the statement: Real and imaginary parts of along the line AB crossing the conducting region is plotted in Figure 8 for FEM and FEM/FEM-Green methods (FEM-BEM is not applied here).  The losses have also been computed. A total of 9.45 mW/m has been obtained with the same mesh. This value should be compared with the value 9.42 mW/m computed with FEM and a circular box with a radius equal to about 100 times the mean size of the device. Those global quantities are quite similar considering the approximations that have been adopted.

Discussion
As a general comment about our method, some simulations have been carried out with a larger inside/outside interface around the active region, i.e., the magnetic circuit in the first example, eventually by including the excitation inside the FEM part. The same tests have been realized with the conducting part in the second example. The results are extremely close to each other in both cases.
However, we notice experimentally that it seems that a regular mesh with a single layer of triangle between the active (FEM) and outside parts gives a better accuracy of the method. Then, some precautions should be taken into account when preparing the mesh with our method. This is a constraint in our method that should be carefully considered at the programming stage.

Conclusions
The FEM/FEM-Green method presented in this paper has been proposed as an alternative to the classical FEM/BEM methods for the treatment of magnetic open boundary problems. The FEM-Green formulation applied to the open region is obtained by using the FEM approximation of fundamental solutions that leads to a discrete equivalent of the second Green's identity. The coupling with the finite element method used for the inner domains appears to be natural due to the same Galerkin approach in both mathematical treatments. Numerical results have pointed out an extremely good agreement between our technique and the standard hybrid FEM/BEM coupling that has been implemented for the sake of comparison. The advantages of the FEM-Green approach are an easier implementation since no particular integration is involved in the computation of the matrix coefficients, the absence of the somewhat awkward normal derivative of the fundamental solution in the formulation, and a more natural coupling with FEM. However, a layer of finite element mesh is required around the inner parts, and it should be carefully taken into account when programming the method. Our numerical scheme can be theoretically extended to 3D problems in the case of scalar potential applications (electrostatics, magnetostatics with scalar magnetic potential), and higher-order (e.g., quadratic) elements. In the latter case, there is no obvious geometrical interpretation of the nodal flux as with linear triangles, but it remains a representative variable. These aspects should be investigated in future work.