Anderson S Nunes

Olivier Chadebec

Patrick Kuo-Peng

Patrick Dular

Gérard Meunier

A Coupling Between the Facet Finite Element and Reluctance Network Methods in 3-D

Keywords: Reluctance network, Finite Element Method (FEM), facet shape functions, Facet FEM (FFEM)

come

I. INTRODUCTION

O PTIMIZING electromagnetic devices can require a large amount of data which, if obtained by means of numerical analysis, can be a time-consuming task. Sometimes, only analyzing large size problems might be time-consuming, mainly when the models have high aspect ratio, for instance, the large electrical machines or actuators with small air gaps.

Many numerical methods [START_REF] Tonti | Why starting from differential equations for computational physics?[END_REF] are available for a wide range of applications. The FEM (Finite Element Method) and the Reluctance Network Method (RNM) are the most widely used for magnetostatic modeling.

The FEM is well known for its flexibility and generality, since the mathematical formulations are solved based on a mesh. Furthermore, the literature available for this method is noticeable, for instance [START_REF] Bastos | Electromagnetic Modeling by Finite Element Methods[END_REF] [START_REF]Magnetic Materials and 3D Finite Element Modeling[END_REF]. However, it can lead to a high number of degrees of freedom and so quite long computation times [START_REF] Turowski | Rapid calculation of stray fields and losses in power transformers[END_REF].

On the other hand, the RNM is one of the most primitive methods for magnetic modeling and its application is based on a reluctance network. This method has remained useful due to its coherent results obtained with low computational effort and computational simulation time [START_REF] Sykulski | Computational magnetics[END_REF]. It has largely been applied to model power transformers [START_REF] Turowski | Rapid calculation of stray fields and losses in power transformers[END_REF], [START_REF] Turowski | Industrial verification of rapid design methods in power transformers[END_REF]- [START_REF] Lopez-Fernandez | 3-D methodology for the heating hazard assessment on transformer covers[END_REF], rotating electrical machines [START_REF] Petrichenko | Contribution à la modélisation et à la conception optimale des turbo-alternateurs de faible puissance[END_REF]- [START_REF] Tang | Analytical modeling of flux-switching machines using variable global reluctance networks[END_REF] and transmission lines [START_REF] Bormann | Reluctance Network Treatment of Skin and Proximity Effects in Multi-Conductor Transmission Lines[END_REF]. In general, the reluctance networks are defined manually, which is usually a hard, long and non-general task.

A methodology for coupling both methods would thus be of importance. Some portions of the studied domain could be described by lumped analytical reluctances, i.e., RNM, and the more complex regions could be automatically discretized and refined as necessary with FEM. This kind of method is known to offer a good trade-off between simulation accuracy and problem size, requiring less time if compared with the fully discretized method [START_REF] Gersem | Field -circuit coupled models in electromagnetic simulation[END_REF].

Manuscript received March 3, 2017. Accepted June 25, 2017. Corresponding author: A. Nunes (email: nunes.anderson@gmail.com).

In [START_REF] Dular | Dual finite element formulations for lumped reluctances coupling[END_REF], a methodology that couples classical nodal/edge 3-D FEM with external reluctance networks is presented. It requires the definition of cuts along the flux tubes to manage the surface jumps of derived scalar potentials as well as the application of a gauge condition using the tree/co-tree technique. Furthermore, in [START_REF] Lehikoinen | A Systematic Approach for Solving Coupled Reluctance Network and Finite Element Models[END_REF], an approach is presented to couple reluctance networks with an FEM magnetic scalar potential formulation, which requires cuts to deal with multiply connected domains.

The conversion of an FEM mesh into an equivalent network has been discussed since [START_REF] Demenko | Reluctance network formed by means of edge element method[END_REF]. It is presented using nodal, edge and facet elements in [START_REF] Demenko | Reluctance network formed by means of edge element method[END_REF]- [START_REF] Demenko | Analogies between finite-difference and finite-element methods for scalar and vector potential formulations in magnetic field calculations[END_REF]. The FFEM (FEM with Whitney facet basis functions) is also applied in [START_REF] Vishnevsky | Boundary integral computation of electric fields in multielectrode galvanic systems using normally continuous elements[END_REF] and [START_REF] Vishnevsky | Modeling magnetization using Whitney facet elements[END_REF] to model static fields and seems to be an appropriate solution for the coupling with the RNM, since its application considering a magnetostatic formulation [START_REF] Vishnevsky | Modeling magnetization using Whitney facet elements[END_REF]- [START_REF] Le-Van | Développement de formulations intégrales de volume en magnétostatique[END_REF] generates an FFEM stiffness matrix that can be interpreted as a magnetic circuit [START_REF] Araujo | Plateforme de conception d'actionneurs: méthodes et outils pour le pré-dimensionnement d'actionneurs[END_REF] [START_REF] Araujo | 3-D hybrid FEM & BEM using whitney facet elements and independent loops[END_REF].

Considering all the complexities of applying the classical 3-D FEM formulations in complex domains and coupling these formulations with external reluctance networks, this paper aims to present an alternative solution. Once the FFEM matrix can be interpreted as a magnetic circuit and easily coupled with the RNM, this entire system, fully compatible with both numerical approaches, i.e., FEM and RNM, can be solved by a single 0-D circuit solver that unifies the magnetic induction free divergence constraint. This network can be solved considering magnetomotive force (mmf ) sources from permanent magnets, coils or from external magnetic circuits. The flux sources are allowed in the external circuits and can also be used to impose magnetic induction fields along the boundaries of the discretized domain.

For this purpose, the FFEM formulation, the way to apply its boundary conditions (BCs) following the circuit formalism and the coupling of the network obtained through the FFEM with the external one are presented. The nodal and independent loop circuit approaches are analyzed and the FFEM formulation is compared with the classical edge A-formulation.

Finally, one simple case is used to illustrate the FFEM matrix and the coupling process. Then, the results of modeling an actuator with the classical FEM A-formulation and with the proposed methodology are compared.

II. MAGNETOSTATIC FACET FEM METHODOLOGY A. Basis equations

The magnetostatic fields are described by the Ampere (1) and Gauss (2) equations, together with the magnetic constitutive relation for linear and isotropic materials (3), i.e., curl

H = J s (1) div B = 0 (2)
B = µH + B r ( 3
)
where H is the total magnetic field, J s is the source current density, B is the magnetic induction, µ is the magnetic permeability and B r is the remnant magnetic induction of permanent magnets.

Total field H can be split into two terms,

H = H s + H r (4)
where reaction field H r is obtained from the reduced magnetic scalar potential ψ,

H r = -grad ψ (5)
and magnetic source field H s , satisfying curl H s = J s [START_REF] Turowski | Industrial verification of rapid design methods in power transformers[END_REF] can be obtained with the Biot-Savart equation, possibly adapted for particular geometries [START_REF] Urankar | Vector potential and magnetic field of current-carying finite arc segment in analytical form, part II: Thin sheet approximation[END_REF]- [START_REF]Vector potential and magnetic field of current-carrying circular finite arc segment in analytical form -Part V. Polygon cross section[END_REF]. The source of field H s can thus be outside the studied domain. If it is partially or completely inside the domain, conformal meshes are not required. Besides that, it is possible to use higher order elements for the integration process, increasing the accuracy of H s . Applying ( 5) in ( 4), one has

H + grad ψ = H s (7)
Applying (3) in [START_REF] Soto | Software for fast interactive three-dimensional modeling of electromagnetic leakage field and magnetic shunts design in Shell type transformers[END_REF], one obtains the strong formulation based on B, i.e.,

1 µ B + gradψ = H s + 1 µ B r (8)

B. Facet shape functions

The shape function spaces are defined using the Whitney forms [START_REF] Whitney | Geometric integration theory[END_REF] [START_REF] Bossavit | Whitney forms: a class of finite elements for threedimensional computations in electromagnetism[END_REF], where W 0 , W 1 , W 2 and W 3 are the nodal, edge, facet and volume function spaces, respectively. The relations between these spaces are given by the Rham's diagram [START_REF] Bossavit | Whitney forms: a class of finite elements for threedimensional computations in electromagnetism[END_REF], The facet shape function is plotted for the face parallel to the XZ-plane on a reference tetrahedral element, where the dimensions are unitary and the referred face has an area of 0.5. Its behavior is depicted in Fig. 1a and its normal component is presented in Fig. 1b.

W 0 grad f ---→ W 1 curl F ---→ W 2 div F ---→ W 3 (9)
As illustrated in Fig. 1b, the flux is equal to 1 for a given face i

W i • n i = 1 S i (10)
and null along the other faces, Fig. 1a. S i is the surface area of face i.

Another important aspect of the facet function is that its divergence in an element is equal to the inverse of the volume V e of this element (for first order elements):

div W i = 1 V e (11)
C. Magnetostatic formulation based on facet shape functions Projecting (8) on a domain Ω (Fig. 2) by the Galerkin method, with the facet shape functions B ∈ W 2 as test functions, leads to ( 12)

Ω 1 µ B • B dΩ + Ω grad ψ • B dΩ - Ω H s • B dΩ- Ω 1 µ B r • B dΩ = 0 (12)
Since the magnetic induction B is given by

B = n f i=1 B i Φ i (13)
for first order elements, the first term of the left side of ( 12) can be rewritten as:

Ω 1 µ B • B dΩ = n f j=1   Ω 1 µ B j • B i dΩ   Φ i ( 14
)
where j and i are the faces indices, n f is the number of faces and Φ i is the magnetic flux through the face i. It is important to notice that ( 14) results in a self-reluctance for i = j and in a mutual reluctance between faces i and j for i = j.

The second term of ( 12) is solved using the following vector identity div

(f A) = grad f • A + f div A (15)
resulting in:

Ω grad ψ • B dΩ = Ω div (ψ B ) dΩ - Ω ψ div B dΩ (16)
Applying the divergence theorem in the second term of ( 16), one obtains

Ω grad ψ • B dΩ = Γ ψ B • n dΓ - Ω ψ div B dΩ (17)
where Γ is the boundary of Ω, Fig. 2.

The surface integral term of ( 17) is null and needs to be analyzed in the following cases (See Fig. 2): Thus, considering that inside each element Ω l or Ω m , the reduced scalar potential ψ l or ψ m are constants since each of the elements represents one vertice in the dual mesh, i.e. one node in the reluctance network, and also defining the flux direction from Ω l to Ω m , [START_REF] Demenko | Representation of windings in the 3d finite element description of electromagnetic converters[END_REF] becomes

Γ l Face i Γ H Γ B
Ω l +Ωm grad ψ • B dΩ = ψ l Ω l div B dΩ l - ψ m Ωm div B dΩ m (19)
The divergence of B being the inverse of the element volume, as given in [START_REF] Demenko | Representation of windings in the 3d finite element description of electromagnetic converters[END_REF], equation [START_REF]Network models of three-dimensional electromagnetic fields[END_REF] becomes

Ω l +Ωm grad ψ • B dΩ = ψ l -ψ m (20)
Finally, ( 12) is rewritten as

n f j=1   Ω 1 µ B j • B i dΩ   Φ i + ψ a -ψ b = Ω H s • B dΩ + Ω 1 µ B r • B dΩ (21)
or in a matrix form:

[ ] [Φ] -[ψ m ] = [ψ 0 ] (22)
where [ ] is the FFEM matrix, physically known as a reluctance matrix, [Φ] is the unknown flux matrix, [ψ m ] is the magnetic potential between each reluctance and [ψ 0 ] contains the mmf branch sources, defined by means of H s and B r .

Another possibility to obtain this term would be to describe the conductors as stranded current loops, as presented in [START_REF] Demenko | Representation of windings in the 3d finite element description of electromagnetic converters[END_REF]. Even though the magnetic induction free divergence (2) has not been used to obtain the matrix system [START_REF] Demenko | Analogies between finite-difference and finite-element methods for scalar and vector potential formulations in magnetic field calculations[END_REF], it will be considered by solving this matrix system as a magnetic circuit, i.e., a 0-D system where its free divergence is imposed by the Kirchhoff's current law here applied to the magnetic flux.

D. Implementation of n • B| Γ B = 0 BCs

The n • B| Γ B = 0 BC, discussed in Section II-C, implies a zero magnetic flux through surface Γ B . From the circuit standpoint, this means that the branches that represent the element surfaces on Γ B are open. So, they can be suppressed from the equations system, i.e., from the magnetic circuit. This process is exemplified in Fig. 6, where the gray reluctances are not taken into account.

E. Implementation of n × H| Γ H = 0 BCs

The n × H| Γ H = 0 BC can be used to apply symmetry planes and define the BC along coupling surfaces. This is achieved by creating one external node P ( Fig. 3) and connecting it with all the elements along Γ H , which consists in imposing a constant magnetic scalar potential on this boundary, [START_REF] Vishnevsky | Boundary integral computation of electric fields in multielectrode galvanic systems using normally continuous elements[END_REF].

n × H| Γ H = -n × grad ψ| Γ H = 0 (23) Ω Γ H P Fig. 3. BC n × H| Γ H = 0
If the node P is not connected to any other external branch, as presented in Fig. 3, it leads to the magnetic induction free divergence along this surface.

F. Coupling with the external network

The node P (Fig. 3) can be used to connect the meshed domain Ω (converted into a reluctance network) to the external magnetic circuit, as depicted in Fig. 4, where Γ c is the connection surface. As Γ c is defined in the same manner as Γ H (Fig. 3 and ( 23)), the n × H| Γ H = 0 BC is also applied on Γ c .

In terms of implementation, this coupling is performed by adding the external network branches/nodes into the incidence matrix that represents the FFEM system, as well as by adding the external sources and reluctances into their respective matrices presented in [START_REF] Demenko | Analogies between finite-difference and finite-element methods for scalar and vector potential formulations in magnetic field calculations[END_REF].

This appears as an easy process, due to the circuit formalism. Thus, with both networks coupled, i.e., Ω + external network, the entire model can be solved as a single 0-D system.

Ω

G. Circuit modeling

The resulting circuit can be solved using both independent loop and nodal approaches [START_REF] Chua | Computer-Aided Analysis of Electronic Circuits: Algorithms and Computational Techniques[END_REF]. Given a tetrahedral mesh with n elements and e f external surfaces, the circuit characteristics can be analyzed through

b = 2n + e f 2 (24) b = i + n -1 (25)
where b, [START_REF] Vishnevsky | Modeling magnetization using Whitney facet elements[END_REF], is number of branches in the dual mesh, i.e., faces in the mesh [START_REF] Gumhold | Tetrahedral mesh compression with the cut-border machine[END_REF]. Thus, thanks to ( 24) and ( 25), the number of independent loops i can be obtained,

i = n + e f 2 + 1 (26)
So, we may remark that, depending on the BCs, the independent loop approach can generate only a few more unknowns, i, than the nodal approach, which generates n unknowns. Nonetheless, the independent loop approach requires an algorithm capable of defining the tree/co-tree, for instance, the Welsh algorithm [START_REF] Bastos | Electromagnetic Modeling by Finite Element Methods[END_REF] or the method proposed in [START_REF] Nguyen | An independent loops search algorithm for solving inductive PEEC large problems[END_REF]. As mentioned in [START_REF] Le-Van | Développement de formulations intégrales de volume en magnétostatique[END_REF], this process can be timeconsuming in larger meshes. Whereas, the nodal approach does not require any other algorithm to find the tree/co-tree, but it requires the inversion of the reluctance matrix, which can be computationally heavy. Furthermore, since this inversion generates a full matrix, it would impact the performance of sparse matrix tools possibly applied.

After evaluating the solving time, the nodal approach has been used to solve the magnetic circuits, thereby corroborating with the conclusion presented in [START_REF] Demenko | Reluctance network formed by means of edge element method[END_REF]. In the future, this analysis will be extended to nonlinear materials, where the circuit approach impacts on the number of iterations of the nonlinear solvers [START_REF] Derbas | A comparison of nodal-and mesh-based magnetic equivalent circuit models[END_REF].

H. Comparison with A-formulation

The A-formulation based on edge elements is very common in magnetostatic modeling and its comparison with the Bformulation used in this work is of importance. For this purpose, the A weak formulation, with test function

A ∈ W 1 , is Ω 1 µ curl A • curl A dΩ- Ω 1 µ Br • curl A dΩ+ Γ (n × H) • A dΓ = Ω J • A dΩ (27)
Through the functional spaces relation curl A = B and ( 1), ( 27) can be rewritten as:

Ω 1 µ curl A • B dΩ - Ω 1 µ Br • B dΩ+ Γ (n × H) • A dΓ = Ω curl H • A dΩ (28)
Expressing A in terms of B (B = curl A), leads to

Ω 1 µ B • B dΩ - Ω 1 µ Br • B dΩ+ Γ (n × H) • A dΓ = Ω curl H • A dΩ (29)
In the surface integral term, one may rewrite

(n × H) • A as -(n × A ) • H that is null because of the BC n × A | Γ B = 0.
So, considering the boundaries conditions for the magnetic induction B, the formulations are equivalent, since

n × A | Γ B = n • B | Γ B = 0.
Regarding the last term of ( 29), it is possible to apply the Green equation of type curlcurl and curl

A = B , Ω curl H • A dΩ = Ω H • B dΩ + Γ (A × n) • H dΩ (30)
where the last term is null and the source field H s is obtained with the Biot-Savart equation, curl H r = 0,

Ω curl H • A dΩ = Ω H s • B dΩ (31)
Both formulations appear to be equivalent, but considering that they are implemented in different ways, it is worthwhile analyzing others aspects: a) considering the number of edges in a regular tetrahedral mesh [START_REF] Gumhold | Tetrahedral mesh compression with the cut-border machine[END_REF] and the number of unknowns per edge of the A-formulation [START_REF] Takahashi | Investigation of effectiveness of edge elements[END_REF], it was found out that the total number of unknowns are similar in both formulations, b) A-formulation requires tree gauge conditions while Bformulation requires the incidence matrix.

III. MODELS AND RESULTS

The models used to verify the proposed methodology follow an order of complexity. The first one is related to a simple 2-D triangular mesh of an air gap without leakage flux. Then a 3-D trivial model is used to compare the convergence of both formulations. The last two are 3-D magnetic devices, where the magnetic core is simplified using the RNM, i.e. it is modeled by a simple magnetic circuit manually defined, and only the air gap regions are meshed and solved with FFEM. The results are compared with those obtained through FFEM considering the complete model.

The geometries and meshes are generated with the Gmsh program [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and the A-formulation is solved using the GetDP [START_REF] Dular | A general environment for the treatment of discrete problems and its application to the finite element method[END_REF] code.

A. Air gap without leakage flux

As presented in Fig. 5, the 2-D air gap model is defined by one region with the same dimensions and without leakage flux. This is coupled with an external network, where the reluctances are neglected ( = 0). For the sake of simplicity, the relative magnetic permeability is considered unitary and the mmf source is arbitrarily defined as 2500A, resulting in a magnetic flux of 2500Wb. The internal network obtained using FFEM is shown in Fig. 6. This network is then coupled to an external one, resulting in a single circuit. The magnetic circuit is solved and the fluxes values are presented in Table I. The flux values and directions are in accordance with the circuit presented in Fig. 6.

B. Convergence comparison

As presented in Section II-H, the B and A-formulations are similar, consequently presenting the same level of convergence. For the sake of comparison, a simple 3-D model composed by two cubes is created and the dimensions are given in Fig. 7. The inner cube is a PM (permanent magnet) with a coercive field of 890000A/m and relative permeability of 1.09978, while the outer one is a surrounding air domain. The convergence is analyzed considering the magnetic energy calculated for a range of different meshes, from 96 up to 14482 elements, and comparing these values with the energy calculated with a very refined mesh, with 524175 elements. The results are presented in Fig. 8. 1XPEHURIHOHPHQWV

&RQYHUJHQFH>@

)(0$)RUPXODWLRQ ))(0%)RUPXODWLRQ

C. Magnetic device

A more realistic 3-D magnetic device is now considered. It is composed of one winding and two magnetic regions separated by an air gap (Fig. 9). The magnetic permeability of the core is µ = 2000µ 0 and the winding is excited with N I = 3000A. The depth of the magnetic parts, air gap regions, and the outer air domain are 20, 50 and 60mm, respectively. One probe line is placed at mid-thickness of the air gap region. As an additional convergence test, the problem is solved considering a coarse and a fine meshes made of 4151 and 14928 tetrahedral elements respectively. For the fine mesh, the magnetic inductions calculated with both formulations are shown in Fig. 10 and Fig. 11. For the coarse and fine meshes, their Z-components are plotted along the probe line, showing a good matching (Fig. 12). B and A-formulations generate practically the same number of unknowns, 14928 and 14921, respectively, for the fine mesh. The magnetic circuit is then simplified as a coupling of a classical reluctance network and the air gap regions (defined in Fig. 9), the latter asking for a more accurate model, here with FFEM, due to the leakage flux. Two cases with different levels of simplifications are considered (Fig. 13), taking account of some core portions or not:

• coupling 1 (Fig. 13a): some core portions are kept in the FFEM model and their connection surfaces with the RNM model are the core cross-sections intersected by the boundary of the reduced FFEM domain; • coupling 2 (Fig. 13b): the core portions are removed from the FFEM model and their connection surfaces are the core outer surfaces used in coupling 1; only the air regions are discretized. In Fig. 13, the connection surfaces Γ c are represented by continuous red lines and the surfaces with n • B| Γ B = 0 BC are represented by the dashed blue lines. The reluctances 1...4 are for the core portions not considered in the FFEM domains and 5 approximates the winding leakage flux path.

The meshes used for these models (Fig. 14) are clones of the fine mesh of the related portions in the complete model.

Both coupled problems are solved and compared in terms of the relative computational effort and accuracy. The number of unknowns, solution time and the magnetic flux through the upper magnetic part are compared in Tables II and III. In the complete model, the flux is obtained integrating the magnetic induction along the cross section of the upper magnetic part, while the flux through 1 is considered for the problems As a local comparison, the Z-component of B along the probe lines (Fig. 13) is shown in Fig. 15 for all models. The relative differences of this solution between the coupled and full models are shown in Fig. 16. 3RVLWLRQDORQJWKHSUREHOLQH>PP@ %]>7@ 3UREHOLQHLQ WKHDLUJDS 3UREHOLQHLQ WKHDLUJDS

))(0)LQH &RXSOLQJ &RXSOLQJ The differences are small in the region of interest (around 0.1%). However, they increase when the analysis point goes towards the n • B| Γ B = 0 BC, which is a natural consequence since the leakage flux is not totally in the Z-direction. Even though a weak field is computed close to those regions, this behavior cannot be desirable and can be minimized by increasing the X dimensions of the air regions, possibly unifying the domains. A more interesting and definitive solution would be to apply the "infinite" finite elements theory [START_REF] Melgoza | Mapped infinite elements for the solution of open electric and magnetic field problems[END_REF] in the X direction on those surfaces.

As expected, the 'coupling 2' problem presents higher differences, but the results are still coherent and acceptable.

As the solution along the meshed domains is obtained with the FFEM, it is possible to post-process the fields. So, the magnetic induction distribution is obtained in those domains and a zoom in the air gap is presented in Fig. 17. Once the plots are in the same color scale, it is possible to see that the Coupling 2 problem is properly modeling the leakage flux in the air gap. In these coupled models, the winding is completely outside the meshed domain. However, it is possible to consider the windings inside or even partially inside of it. This last possibility is achievable because the SF is calculated independently, in this case using Biot-Savart, and thanks to the circuit coupling. It allows reducing the studied domain even when it does not have symmetry, needing only the connection surfaces. This is an important advantage over the classical formulations.

In order to illustrate this possibility, the model presented in Fig. 9 is modified as shown Fig. 18, where the windings are arranged vertically. The dimensions are the same as the model shown in Fig. 9, except the height of the windings which is 22mm.

X Z

Fig. 18. Geometry -Vertical windings

Because the source field cannot be neglected on the connection surfaces Γ c used in the cut of the upper magnetic part (Fig. 17a), an mmf distribution can appear in the elements adjacent to these surfaces. Consequently, grad ψ| Γ H = 0, which does not respect the assumption defined in Section II-D.

This model is thus simplified as shown in Fig. 19, where the half of the windings are inside the meshed domains and the others halves of the Ampere-turns are considered as mmf sources in the external circuit. The source field is obtained along the entire domain and considering both parts of the windings. This problem is solved with the proposed coupling methodology and the magnetic induction distribution is presented in Fig. 20. Then, it is compared with the results obtained with the FEM A-formulation solved in the complete model, showing a good correlation in the regions of interest (Fig. 21). The number of unknowns was reduced from 19982 (FEM A-formulation) to 13717 (FFEM+RNM).

IV. CONCLUSION

A methodology capable of converting an FEM mesh into a reluctance network, allowing its easy coupling with the RNM, has been presented. As the circuit works as a gauge condition and the formulations are equivalent, it has been possible to keep the accuracy obtained with the A-formulation, but without worrying about the cuts needed to deal with multiply connected domains.

The FFEM analyzed regions have been discretized with a mesh fine enough to ensure the accuracy level and the other outer regions have been simplified using the RNM. Even though only one type of flux tube was used in this work, it allows to use all the range of flux tubes available for this method [START_REF] Nunes | Air Core Reactor Analysis Based on RNM Method[END_REF]- [START_REF] Roters | Electromagnetic devices[END_REF].

Even if the representation of an FEM mesh in terms of a reluctance network has already been shown before, [START_REF] Demenko | Reluctance network formed by means of edge element method[END_REF] for instance, a didactic case (Fig. 6) has been presented in order to show the coupling of this network with an external one, contributing to a different and simpler way to understand the numerical analysis. Additionally, two real cases with different windings arrangements, i.e. horizontal and vertical, were studied. In the last one, only part of the windings was considered inside the studied domain. Such capability allows a wide range of new possibilities to apply this coupling.

With the proposed methodology, the modeling domain can be drastically reduced, meshing only the regions with complex magnetic flux behavior. Two levels of simplifications have been presented, reducing the number of unknowns to 35% and the solving time to 7% in comparison with the complete model, keeping the quality of the results along the region of interest.

This methodology seems to provide more flexibility when dealing with the trade-off simulation time vs. accuracy, mainly for optimizations processes and even when analyzing large models.

It can easily be applied to other static physics with electric circuit analogy, like heat and current conduction [START_REF] Vishnevsky | Boundary integral computation of electric fields in multielectrode galvanic systems using normally continuous elements[END_REF].

Fig. 1 .

 1 Fig. 1. Facet shape functions obtained for tetrahedral element

Fig. 2 .

 2 Fig. 2. The complete domain Ω with its outer surface Γ = Γ H ∪ Γ B and two adjacent elements Ω l and Ωm, with their faces Γ l and Γm.

Fig. 4 .

 4 Fig. 4. Connection between the networks

Fig. 5 .

 5 Fig. 5. Connection between the networks

35 Fig. 6 .

 356 Fig. 6. Meshed domain coupled with the external network

15 20 15

 1520 All dimensions are in mm.

PMFig. 7 .

 7 Fig. 7. Geometry for the convergence test

Fig. 8 .

 8 Fig. 8. Convergence comparison

Fig. 9 .

 9 Fig. 9. Geometry -horizontal winding

Fig. 10 .Fig. 12 .

 1012 Fig. 10. Magnetic induction distribution -FFEM (B-formulation)

Fig. 13 .

 13 Fig. 13. Coupling between FFEM and RNM

Fig. 14 .

 14 Fig. 14. Meshes of the air gap regions (with clipping)

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Bz obtained with the FFEM+RNM coupling

Fig. 17 .

 17 Fig. 17. Magnetic induction distribution along the discretized domains -Zoom in the air gap

3 mmf 2 mmf 1 Fig. 19

 32119 Fig. 19. Coupling FFEM+RNM

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Magnetic induction distribution

TABLE II COMPUTATIONAL

 II EFFORTS OF RNM-FFEM COUPLED MODELS COMPARED TO THE FULL FFEM MODEL

	Method	Number of unknowns	Number of unknowns [%] time [%] Solution
	FFEM (Fine)	14928	Ref.	Ref.
	Coupling 1	8257	55.31	17.09
	Coupling 2	5295	35.47	6.96
		TABLE III
	ACCURACY OF RNM-FFEM COUPLED MODELS COMPARED TO THE FULL
		FFEM MODEL
	Method	Magnetic flux [Wb] Dif. [%] × 10 -4
	FFEM (Fine)	2.5310	Ref.
	Coupling 1	2.520	0.553
	Coupling 2	2.515	0.632

• External surface Γ B : imposes that field B only has tangential components on this surface; • Internal faces of Ω: in accordance with the facet functions definitions, B is null everywhere except inside the two elements that share the face i; • Γ l ∪ Γ m : This second case allows to limit Ω to the sum of Ω l to Ω m , Fig.2and (18). However,B • n is null along Γ l ∪ Γ m , except along Γ l ∩ Γ m ; • Γ l ∩ Γ m :The shape functions B l and B m have opposite directions along these faces, canceling each other out.Ω l +Ωm gradψ • B dΩ = Γ l +Γm ψ B • n Γ dΓ-Ω l +Ωm ψ div B dΩ(18)nΩ l Ω m Γ m Γ Ω

ACKNOWLEDGMENT This work is supported by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico -Brazil), project 200215/2016-7. Anderson Nunes acknowledges Engineering Simulation and Scientific Software Ltda. (ESSS) for the help during this work.


Anderson Santos Nunes

Patrick Dular

Olivier Chadebec

Patrick Kuo-Peng

Anderson S Nunes

Subproblems Applied to a 3-D Magnetostatic Facet FEM Formulation

Keywords: Reluctance Network Method (RNM), Facet Finite Element Method (FFEM), Magnetic Devices, Three-Dimensional, Magnetostatics

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The Facet Finite Element Method (FFEM) has been presented as an alternative to reduce the computational time and resources when modeling magnetostatic fields in electromagnetic devices [START_REF] Nunes | Magnetic devices analysis by Face FEM coupled with standard reluctance network[END_REF]- [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF], because it allows an easy coupling with external magnetic networks, i.e., the Reluctance Network Method (RNM), and does not require any additional attention to deal with multiply-connected domains as in the classical magnetic scalar potential formulation.

However, the FFEM magnetostatic formulation does not take into account an electrical current source, requiring some additional computations when the discretized domains contain windings. In [START_REF] Demenko | Network models of three-dimensional electromagnetic fields[END_REF] and [START_REF] Demenko | Analogies between finite-difference and finite-element methods for scalar and vector potential formulations in magnetic field calculations[END_REF] the sources are computed considering the windings as current loops. A methodology using the electric vector potential T0 is presented in [START_REF] Wojciechowski | Description of the windings of the electromagnetic energy converters using the modified T0 method[END_REF].

Considering all the complexities of calculating these sources, a possible solution could be to take advantage of subproblem method (SPM), which has been applied to model electromagnetic problems in a progressive way [START_REF] Dular | Subdomain perturbation finite element method for skin and proximity effects in inductors[END_REF], allowing to benefit from a previous solution, instead of starting a new complete one for any physical variation, e.g., material properties or geometric.

This methodology has been applied to model nondestructive testing [START_REF] Badics | An effective 3-D finite element scheme for computing electromagnetic field distortions due to defects in eddycurrent nondestructive evaluation[END_REF], to take into account the effects of air gaps [START_REF] Dular | Perturbation finite element method for magnetic model refinement of air gaps and leakage fluxes[END_REF] and in general magnetodynamics applications [START_REF] Dular | Subdomain finite element method for efficiently considering strong skin and proximity effects[END_REF]- [START_REF] Dular | Progressive current source models in magnetic vector potential finite-element formulations[END_REF] applying the A or the A-V formulation. It was also used to model MEMS through an electrostatic formulation, where an unperturbed electric field is obtained and then a conductive domain is inserted as a perturbation [START_REF] Boutaayamou | A perturbation method for the 3D finite element modeling of electrostatically driven MEMS[END_REF].

Thus, considering that i) normally the authors start from an unperturbed solution and then apply all the perturbations necessary to take into account the physical phenomena and/or geometric characteristics in interest, ii) the source field needed in the FFEM magnetostatic formulation can be obtained with the Biot-Savart equation which can be interpreted as an unper-Manuscript received December 1, 2012; revised August 26, 2015. Corresponding author: A. Nunes (email: nunes.anderson@gmail.com).

turbed solution, and iii) all these SPM magnetic applications above cited are using A based formulations, this paper aims to present the SPM applied to the FFEM B-formulation, reducing the effort required to obtain the source field solution.

This is accomplished by considering the Biot-Savart solution as a first subproblem, i.e., the unperturbed one. Then, the active parts (µ r > 1) are added into the domain and a reaction field solution is obtained with the FFEM considering the Biot-Savart source field only inside these ferromagnetic regions. The total solution is the sum of both problems (unperturbed + perturbation). This is the key point of this work, since it allows reducing the source field calculation to the active parts only, decreasing the computational time and effort.

Furthermore, it is presented how to correct the boundary condition difference between the FFEM and Biot-Savart, which allows a direct superposition of the results.

In order to validate the obtained results, two test cases are arranged. The first one is a simple limit case built based on a 2-D mesh, which was extruded creating the 3-D domain. The second test case is a representative 3-D magnetic device that is modeled using the proposed methodology. Finally, FFEM + RNM coupling [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF] and the SPM are applied in this model at the same time, benefiting from the advantages of both methodologies.

Considering that this work deals with magnetostatic fields, the results obtained are compared with those calculated using the classical FEM A-formulation or with the FFEM Bformulation.

II. MAGNETOSTATIC FFEM FORMULATION AND ITS BOUNDARY CONDITIONS

The magnetostatic problem is defined along a bounded Euclidean domain Ω with its boundary Γ = Γ H ∪Γ B (Fig. 1a). Ω a can represent a perfect conductor (σ → ∞) Ω pc , i.e., a domain where internal magnetic field is zero and results in a zero skin depth [START_REF] Dular | Subproblem finite-element refinement of inductors from wire to static and dynamic volume models[END_REF], [START_REF] Dular | Progressive eddy current modeling via a finite element subproblem method[END_REF], or a magnetic material with finite permeability Ω m . Ω s is a stranded conductor in which the current density J s is imposed. The interface between two different media in Ω is denoted by γ (Fig. 1b). The Ampère and Gauss equations, together with the magnetic constitutive relation for linear and isotropic materials, are expressed as

Γ H Γ B Ω a Ω s (a) Domain Ω defined by its boundary Γ = Γ H ∪ Γ B . Ω 1 Ω 2 γ n (b) Interface γ between two me- dia.
curl H = J s (1) div B = 0 (2)
B = µH + B r (3)
or

H = νB + H c ( 4
)
where H is the magnetic field, B is the flux density, µ is the magnetic permeability, ν is the magnetic reluctivity, B r is the remnant flux density and H c is the coercive field. In the SPM, the B r or H c can also be used as volume sources, expressing changes in the permeability. The boundary conditions are

n × H| Γ H = 0 (5) n • B| Γ B = 0 ( 6
)
where n is the unit normal vector. The interface conditions between two media in Ω, γ in Fig. 1b, are

n × (H 1 -H 2 ) | γ H = J f (7) n • (B 1 -B 2 ) | γ B = 0 (8)
which can be expressed through the following notation

[n × H] γ H = J f (9) [n • B] γ B = 0 (10)
where J f is the surface current density. γ H and γ B are interfaces between two media considering the interface conditions for H and B fields.

A. FFEM formulation

As previously cited, this work's goal is to apply the SPM in a FFEM magnetostatic formulation, it is briefly summarized here. The details of this formulation can be found in [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF].

It is established writing the total magnetic field H in terms of source and reaction fields, H s and H r respectively,

H = H s + H r (11)
with curl H s = J s [START_REF] Dular | Magnetic model refinement via a perturbation finite element method from 1D to 3D[END_REF] and

H r = -grad ψ ( 13
)
where ψ is the reduced (single valued) magnetic scalar potential. The source field H s [START_REF] Dular | Magnetic model refinement via a perturbation finite element method from 1D to 3D[END_REF] is calculated by applying the Biot-Savart equation, which can be a time consuming task. So, applying [START_REF] Dular | Subdomain perturbation finite-element method for skin and proximity effects[END_REF] and ( 13) in (3) leads to the strong formulation of the magnetic constitutive relation based on B,

1 µ B + grad ψ = H s + 1 µ B r (14)
The weak form of ( 14) is

n f i=1   Ω 1 µ B j • B i dΩ   Φ i + ψ a -ψ b = Ω H s • B dΩ + Ω 1 µ B r • B dΩ (15)
with its boundary condition

n • B| Γ B = 0 (16)
where B ∈ W 2 are the facet shape functions [START_REF] Raviart | A mixed finite element method for 2-nd order elliptic problems[END_REF], [START_REF] Bossavit | A rationale for "edge-elements" in 3-D fields computations[END_REF]. ψ a and ψ b are the magnetic potentials at two adjacent elements a and b. These basis functions have some important properties, for instance, its normal component is constant along the first order element faces. Furthermore, its flux is equal to 1 for a given face i

B i • n i = 1 S i (17)
and zero along the other faces. S i is the surface area of face i.

Another important aspect is that its divergence in an element is equal to the inverse of the volume V e of this element (for first order elements):

div B i = 1 V e (18)
Thus, the final matrix system

[ ] [Φ] -[ψ m ] = [ψ 0 ] + [Br] (19)
is composed by the terms of ( 15): the FFEM stiffness matrix, [ ], physically known as a reluctance matrix; the unknown magnetic flux vector, [Φ]; the averaged magnetic potential between each reluctance, [ψ m ]; the vector with all the magnetomotive force (mmf ) sources, [ψ 0 ], obtained considering the magnetic field H s calculated using Biot-Savart; and the vector which takes into account any remnant induction, [Br] . Since the matrix system is stated without taking into account (2), the zero divergence of B is enforced by ( 19) as a circuit system, where the Kirchhoff's current law is imposed. It suggests that the solution can be obtained by the use of a 0D circuit solver, as presented in [START_REF] Demenko | Network models of three-dimensional electromagnetic fields[END_REF], [START_REF] Demenko | Network equivalents of nodal and edge elements in electromagnetics[END_REF], [START_REF] Demenko | Magneto-electric network models in electromagnetism[END_REF].

As aforementioned, this formulation is detailed in [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF].

B. The n × H| Γ H = 0 boundary condition

Going a step forward in the boundary conditions analysis, besides the n • B| Γ B = 0 (16), it is also possible to apply n × H| Γ H = 0. Considering that ( 19) is solved in the light of the circuit formalism, n × H| Γ H = 0 is achieved by connecting all the elements along Γ H with an external node P (Fig. 2), which consists in imposing a constant magnetic scalar potential on Γ H , as

n × H| Γ H = -n × grad ψ| Γ H = 0 (20) Γ H P Fig. 2. Boundary condition n × H| Γ H = 0.
If node P is not connected to any other external branch, like presented in (Fig. 2), it implies a zero total magnetic flux along this surface.

This condition was discussed in our previous paper [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF], but this connection is also used to impose one external flux density in the model, i.e., n

• B| Γ B . C. Normal flux density ( n • B| Γ B ) source
Considering the FFEM circuit formalism, this source can be implemented defining flux sources along each face i along Γ B , i.e., one flux source at each branch connecting the face i to the node P shown in Fig. 2. The flux values are given by

Φ i = Γi n • B i dΓ i (i = 1..., n f Γ ) (21)
where n f Γ is the number of element faces along Γ and B i is the flux density field along each face. This is not a time-consuming process, since the integration given in ( 21) is based on the Gauss points of those 2-D elements along Γ.

III. DEFINITION OF THE SUBPROBLEMS

Once the key concepts necessary to define the SPM in the light of the FFEM (circuit formalism) are above presented, it is possible to describe the SPM chain.

It is based on a sequence of subproblems, considering the idea of the source and reaction fields, but only calculating the source field along specific parts of the domain.

In this case, a sequence of two problems, p and q, is presented. In problem p an unperturbed field is obtained with Biot-Savart. In problem q, the source field is obtained, also with Biot-Savart, and the reaction field is solved with the FFEM. Then, both results are simply superposed, resulting in the total solution (p + q). Depending on the case, the first and consequently the sequence of the problems can vary. Here problem p is always the Biot-Savart solution and problem q can differ depending on the application.

A. Inductor alone

The Biot-Savart field calculation in subproblem p is obtained in Ω considering only the conductors, Ω s,p , where Q ∈ Ω p is the calculation point, P ∈ Ω s,p is the integration point and r P Q is the position vector, with the current density J s (Fig. 3). For 3-D cases, the source field can be obtained through

H(Q) r P Q Q J s dv, Idl, Jds P Ω sv Ω ss Ω sl
H (Q) = 1 4π Ωsv J (P ) × r P Q |r P Q | 3 dΩ sv (22)
Moreover, for some specific cases where the aspect ratio between the entire domain and conductor cross-section dimensions is high, Ω s,p can be simplified by its average line, not requiring a 3-D but an 1-D mesh which could represent computational advantages. In this case, the source field is obtained using

H (Q) = I 4π Ω sl dl (P ) × r P Q |r P Q | 3 dΩ sl (23)
In case of 2-D symmetries, ( 22) can be re-written as

H (Q) = J 2π Ωss ds (P ) × r P Q |r P Q | 2 dΩ ss (24)
The integration process of ( 22), ( 23) or [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF] along Ω s,p can be done based on a non-conforming mesh in Ω, allowing the utilization of higher order isoparametric elements, for instance, increasing the accuracy of the source field, mainly when Ω s,p contains curved regions.

The source field can be calculated normalizing the current in the conductors as 1A, for instance, and multiplying the field by the real currents afterwards. Also, as the conductors are immersed in a region with µ r = 1, it is possible to explore symmetries, which can be lost in the complete problem [START_REF] Dular | Progressive inductor modeling via a finite element subproblem method[END_REF].

B. Correction of the boundary condition not imposed by the Biot-Savart equation

The Biot-Savart solution, Section III-A, does not impose the boundary condition n • B| Γ B = 0 along the boundaries as the FFEM does. So, a simple superposition of results, i.e., problem p (Biot-savart) + problem q (FFEM), might cause some differences in the total solution. This problem is also verified when using the classical FEM with A formulation, because it also imposes n • B| Γ B = 0, meaning that this correction has to be applied there, as well.

It can be solved canceling the component n • B p along the outer surface of Ω, Γ, using the n

• B| Γ B surface sources presented in Section II-C, with n • B q | Γ = -n • B p | Γ .
As mentioned in Section II-C, it is a light process since it depends on a 2-D integration along those element faces which comprise the surface. In addition, this step can be included in a subsequent problem, q.

C. Perfect conductor boundary condition

Other application of the methodology explained in Section III-B, i.e., the implementation of n • B q | Γ = -n • B p | Γ , is the limit case of a perfect conductor. Furthermore, it is a limit case that, in the light of the FFEM, applies a n • B| Γ field along the boundary as described in Section II-C. So, it appears as an interesting case to be analyzed.

The internal magnetic field of a perfect conductor, denoted by Ω pc ∈ Ω and its boundary γ pc , is zero and results in a zero skin depth. This behavior is treated removing Ω pc from Ω and fixing a zero total magnetic flux along its boundaries, γ pc . So, assuming that there is no discontinuity along γ pc in the total solution, as in [START_REF] Dular | Subdomain finite element method for efficiently considering strong skin and proximity effects[END_REF], it is possible to perturb the problem p, solved with Biot-Savart (Section III-A), defining a problem q as

[n • B q ] γpc,q = B f,q (25)
where B f is a surface field that defines the flux density discontinuity. But, as stated in [START_REF] Dular | Subdomain perturbation finite-element method for skin and proximity effects[END_REF], the total solution has no B f ,

B f = B f,p + B f,q = 0 so, B f,q = -B f,p (26)
consequently,

B f,q = -[n • B p ] γpc,p (27)
As the internal field of Ω pc is null, this discontinuity is treated as a boundary condition,

n • B q | γpc = -n • B p | γpc (28)
As the field B p is only necessary along γ pc , the Biot-Savart field is only solved along this boundary in problem p, reducing the computational time of the entire problem. The total solution, p + q, can also be interpreted as the Dirichlet boundary condition, i.e., n • B| γpc = 0, which imposes a zero normal component on the fields, offering a solution for further comparisons.

D. Volume source

As previously mentioned, the aim of this work is to take advantage of the SPM in order to reduce the domain where the Biot-Savart field is obtained and this is achieved by considering volume sources only along the active parts of this domain.

So, given a first problem p, with its constitutive relation and Maxwell equations,

H p = 1 µ p B p (29) rot H p = J (30) div B p = 0 (31)
it is possible to perturb this problem changing the magnetic permeability of a certain region, from µ p to µ q , where q is a subsequent problem, resulting in a new constitutive relation for the total problem p + q,

H p + H q = 1 µ q (B p + B q ) (32)
The difference between the solutions (29) and (32) gives

H q = 1 µ q - 1 µ p B p + 1 µ q B q (33)
and suggests that problem q depends on the solution of problem p only along those perturbed regions, from µ p to µ q , i.e., the active regions. Moreover, the total solution should satisfy the Ampère equation,

rot (H p + H q ) = J ( 34
)
where it is possible to observe that H q is curl-free. Applying the Gauss equation, div (B p + B q ) = 0 (35) one can see that B q is divergence-free. This means that the field H q has the same behavior of the field generated by a permanent magnet, H c in (4), and allows to take it into account in the B r term of the FFEM magnetostatic formulation given in [START_REF] Dular | Progressive current source models in magnetic vector potential finite-element formulations[END_REF], thanks to (3) and ( 4). Finally, re-writing (33) and defining the total field of the second problem in terms of the source and reaction fields, H q = H qs + H qr [START_REF] Dular | Subdomain perturbation finite-element method for skin and proximity effects[END_REF], one obtains an equation with the same shape of ( 14), which consequently can be solved through the FFEM,

1 µ q B q + grad ψ q = H qs - 1 µ q - 1 µ p µ 0 H p (36)
and contains the volume source term,

VS = - 1 µ q - 1 µ p µ 0 H p (37)
The field H qs is zero because this problem does not consider the source field but only the volume sources.

As the perturbation regards the change from µ p to µ q , there would have no problem to consider non-linear magnetic materials, once µ q can be considered as element based.

If the model contains external reluctances with mmf sources, simplifying possible windings, they are not taken into account in the discretized domain. These mmf sources need to be corrected using the factor given in (37).

Reference [START_REF] Dular | Finite Element Magnetic Models via a Coupling of Subproblems of Lower Dimensions[END_REF] presents a way to correct problem q for higher values of µ q , i.e., H p approaches to -H q . This correction is important in the case of projections of curl A [START_REF] Geuzaine | A galerkin projection method for mixed finite elements[END_REF] between different meshes and that is not the case here because the H p field is calculated directly along the mesh used to solve problem q.

IV. APPLICATIONS

In order to validate the proposed methodology, some test cases are presented following the order of complexity. First, the correction for the Biot-Savart boundary condition is solved together with a perfect conductor simplification on a simple geometry. Then, one 3-D magnetic device is completely solved using this proposed methodology. Finally, this device is simplified using the FFEM+RNM coupling, as proposed in [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF], and solved considering the subproblem technique.

The geometries and meshes are generated with the Gmsh program [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and the A-formulation is solved using the program GetDP, which abbreviation means a "General Environment for the Treatment of Discrete Problems" [START_REF] Dular | A general environment for the treatment of discrete problems and its application to the finite element method[END_REF]- [START_REF] Geuzaine | GetDP: a general finite-element solver for the de Rham complex[END_REF].

A. Correction of the boundary condition not imposed by Biot-Savart

These preliminary test cases are based on a 3-D geometry obtained by the extrusion of a 2-D model, composed of two conductors and a magnetic core, as depicted in Fig. 4a. The top view of the 3-D mesh generated is shown in Fig. 4b. One horizontal probe line placed in the center of the model is used to compare the results locally. In order to show the correction effect solely, the normal component of the magnetic field is obtained along the external surface with Biot-Savart, (Fig. 5), and the relative permeability of the magnetic core is set as unitary. The problem is solved considering the correction presented in Section III-B and the correction field obtained is shown in Fig. 6. The color scale is set as logarithmic in order to show a wider range of field magnitudes. As defined in Section III-B, the direction of the fields H, Fig. 5, and B, Fig. 6, have opposite directions along the outer surface Γ, where the correction is applied.

B. Perfect conductor boundary condition

The magnetic core is now considered as a perfect conductor, Ω pc , and the correction presented in Section III-C is applied along its boundary. The Biot-Savart correction is also considered.

The correction field is shown in Fig. 7 and the comparison with the classical FEM along the probe line is presented in Fig. 8, where B p (BS) is the source field obtained with Biot-Savart, B q is the correction field and B q + B p is the sum of both quantities. One can observe that the field obtained from Biot-Savart, B p (BS), is a continuous solution and when it is summed up to the one obtained with a first order FFEM, the total solution, i.e., B q + B p becomes more continuous. Since this is not the case for the FFEM solution, this fact can be interpreted as a gain of accuracy of the proposed method.

C. Volume sources

Now, a more realistic magnetic device is simulated. As this work is a sequence of [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF], to make the computational time comparison based on the conventional FFEM easier and clearer, the magnetic device presented in [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF] is utilized here as well. Furthermore, its vertical windings close to the air gaps generate higher magnetic fields along these regions, consequently, increasing the problem p contribution in the total solution, which is an important aspect due to the different behaviors of both solutions.

It is composed by two windings fed in opposite directions by a current density of 2.73 10 7 A/m 2 and two magnetic regions, with µ = 2000µ 0 , separated by an air gap, as shown in Fig. A tetrahedral mesh of the active parts and windings is shown in Fig. 11. This magnetic device is completely solved using the volume source described in Section III-D, and the Biot-Savart correction defined in Section III-B.

The field H p presented in Fig. 12 was obtained by Biot-Savart only along the active parts, reducing the computational time to 34% compared to the computational time of the classical solution, where Biot-Savart field is needed along the entire domain. Then, the problem q solution is obtained solving (36) with FFEM [START_REF] Bossavit | A rationale for "edge-elements" in 3-D fields computations[END_REF], considering the source field H p given in Fig. 12 as a volume source, (37).

Even if not necessary to calculate the p problem solution along the entire domain, but only along the interest regions, it is shown in Fig. 13a as a post-processing task. Therefore, the total (p + q) solution of this problem is represented by the sum of these fields (Fig. 13a + 13b), where Fig. 13b shows the solution of problem q. A comparison between the total solution, B p + B q , and the results obtained with the classical FFEM along the probe line is given in Fig. 14. It is important to mention that the problem p consists of calculating the Biot-Savart field along this line, as a post-processing operation.

D. Volume sources + Coupling (FFEM+RNM)

The presented in Fig. 10 is simplified as shown in Fig. 15, where half of the windings are inside the meshed domains and the other halves of the Ampere-turns are considered as mmf sources in the external circuit. The source field, H p , is obtained considering the Biot-Savart integration along the complete windings only along the active parts of the meshed domain, reducing the computational time needed for the source field calculation to about 12% of the computational time of the classical method.

Finally, the volume source (37) is used as a source term in the FFEM system and its correction is also applied in the external mmf sources. So, this coupled problem is solved using also the subproblem technique and the flux density distribution obtained is presented in Fig. 16b. The problem p solution is shown in Fig. 16a, only with the visualization purpose as already mentioned, composing the total solution (Fig. 16a + Fig. 16b). Then, this coupled model is also solved using the classical coupling [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF] and the comparison along the probe line is presented in Fig. 17.

3RVLWLRQDORQJWKHSUREHOLQH>PP@ %]>7@

))(0FRXSOHG5HI B p + B q FRXSOHG Fig. 17. Field Bz along the probe line.

E. Time comparison of the source field calculation

In reference [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF], the source field is computed along the entire domains, but in this work it is calculated only along the ferromagnetic parts. As the test cases are the same for both works, it is possible to compare the simulation time necessary to obtain the source field in all those test cases. The comparison is presented in Table I, where the percentage values refer to the computational time necessary to obtain the source field along the complete model (indicated with "Ref." in Table I) presented in [START_REF] Nunes | A coupling between the facet finite element and reluctance network methods in 3-D[END_REF]. Notwithstanding the coupling FFEM+RNM can save computational effort and time, the results obtained in this work make evident how the subproblem technique can represent an important speedup.

V. CONCLUSION

The proposed methodology allowed to reduce the computational time to obtain the source field required in the Facet Finite Element Method by applying subproblem technique. It was based on the application of the volume source together with a correction for the boundary condition not imposed by Biot-Savart but imposed by the FEM technique, allowing a direct superposition of those two problems, p and q. Computationally speaking, it was possible to find that this correction is not a time consuming task. However, it can represent a significant impact on the accuracy of solution, mainly near the outer surfaces of the domain.

The solutions were presented in the order of complexity, showing the different approaches independently, which can be helpful to apply in others situations. The results obtained for one simple and other complex cases were compared with the results calculated through the classical FEM (A-formulation) and FFEM, presenting an acceptable accuracy.

For the real case studied, it was possible to reduce the Biot-Savart computational time to 12% if compared with the classical method. As this factor is directly dependent on the volume ratio between the magnetic and non-magnetic regions in the domain, it can vary considerably. Nonetheless, normally the 3-D cases demands a considerable large air regions due to the imposition of the n • B| Γ B = 0 boundary condition, this factor tends to be small in the majority of practical cases.

Fig. 1 .

 1 Fig. 1. Domain Ω and the interface γ.

Fig. 3 .

 3 Fig. 3. Integration domain for the Biot-Savart equation.

  2-D geometry (extruded to build a 3-D mesh). (b) Top view of the 3-D mesh.

Fig. 4 .

 4 Fig. 4. Geometry and mesh of the extruded domain.

Fig. 5 .

 5 Fig. 5. Magnetic field obtained with Biot-Savart along the Γ.

Fig. 6 .

 6 Fig. 6. Flux density distribution along Ω.

Fig. 7 .

 7 Fig. 7. Flux density distribution -problem q.

Fig. 8 .

 8 Fig. 8. Field By along the probe line.

Fig. 9 Fig. 9 .

 99 Fig.9shows a close view of By.

10 .Fig. 10 .

 1010 Fig. 10. Geometry of the 3-D device.

Fig. 11 .

 11 Fig. 11. 3-D mesh of the device.

Fig. 12 .

 12 Fig. 12. Source field Hp obtained only along the active parts.

  (a) Magnetic field Hp obtained with BS. Bp = µ 0 Hp. (b) Flux density Bq obtained with FFEM.

Fig. 13 .

 13 Fig.[START_REF] Dular | Progressive inductor modeling via a finite element subproblem method[END_REF]. Field maps representing the solutions of the problems p and q for the complete model.

Fig. 14 .

 14 Fig. 14. Field Bz along the air region.

3 mmf 2 mmf 1 Fig. 15 .

 32115 Fig. 15. Coupling between FFEM and RNM.

  (a) Magnetic field Hp obtained with Biot-Savart. Bp = µ 0 Hp. (b) Flux density Bq obtained with FFEM.

Fig. 16 .

 16 Fig.[START_REF] Boutaayamou | A perturbation method for the 3D finite element modeling of electrostatically driven MEMS[END_REF]. Field maps representing the solutions of the problems p and q for the reduced model.

TABLE I COMPARISON

 I OF THE SIMULATION TIME IN RELATION WITH THE COMPLETE MODEL SOLVED WITHOUT THE SPM.

	Model	Without the SPM	With the SPM
		(From reference [3]) (Proposed method)
	Complete	Ref.	33.6%
	Coupled (FFEM+RNM)	52.2%	12.1%

ACKNOWLEDGMENT

This work is supported by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico -Brazil), project 200215/2016-7.