An efficient numerical scheme for the FE-approximation of magnetic stray fields in infinite domains

In this contribution we propose an efficient and simple finite-element procedure for the approximation of open boundary problems for applications in magnetostatics. In these problems, the interaction of the solid with external space plays a crucial role because of the magnetic stray fields that arise. For this purpose, the infinite region under consideration is approximated by a sufficiently large domain. This region is then divided into a so-called interior domain and an exterior domain. As an essential prerequisite, we assume linear behavior of the (large) exterior domain. The latter is then reduced to the degrees of freedom of the connecting line (2D)/connecting surface (3D) of both domains via static condensation. The proposed finite element scheme can be seen as an alternative to established methods for infinite domains. These methods often require semi-analytical solutions to describe the behavior in the exterior domain, which can be difficult to obtain if heterogeneous structures are present. The proposed finite element procedure is not subject to any restrictions with regard to the topology of the exterior space. After a general introduction of the numerical scheme, we apply the method to problems of magnetostatics with nonlinear behavior in the interior domain.


Introduction
A variety of engineering problems require the solution of open boundary problems, e.g. the simulation of acoustic phenomena or seismic waves. In magnetostatics of solids (magnets), the demagnetization fields inside a magnetic body and the stray fields outside the magnetic solid are directly affected by the magnets. For simple geometries, the stray field, which is proportional to the magnetization, can be captured by demagnetization coefficients. However, more complex geometries usually require a more sophisticated analysis of the considered fields over the whole domain. In B Jörg Schröder j.schroeder@uni-due.de Maximilian Reichel maximilian.reichel@uni-due.de Carolin Birk carolin.birk@uni-due.de 1 these cases, the balance equation to be solved is the Gaussian law, which is defined over the magnetic solid and the surrounding external space. Thus, the analysis of an open boundary problem is necessary. Since the stray field vanishes at a certain distance from the magnet, an approximate analysis can be restricted to a sufficiently large outer spatial region. Defining the boundary of the finite but sufficiently large exterior space is initially an intuitive process. However, a sufficiently precise choice can usually be made based on experience. In a complex application, the homogenization of magneto-elastic problems in the presence of large deformations, was solved in [1,2] by Keip and Rambausek by defining a macroscopic boundary value problem considering a large but finite exterior domain. An overview of finite element open boundary techniques can also be found in [3]. Here, the authors divide the various techniques into 12 subgroups, such as truncation of outer boundaries, ballooning, infinite elements, infinitesimal scaling, absorbing boundary conditions, and others. A rule of thumb for the truncation radius is given in [3]. However, in our case the analysis of the Gaussian law requires a boundary condition at infinity to ensure the well-posedness of the problem, see [4]. A truncation of the outer boundary based on a finite element approach requires the replacement of the condition at infinity by a boundary condition on a finite artificial surface. Estimates of the error due to the boundary values at the finite boundary as well as higher-order boundary conditions to reduce the error are derived in [4] for several applications. A finite-difference time-domain free-space simulation method for solving unbounded electromagnetic problems based on an absorbing layer that absorbs electromagnetic waves without reflection has been proposed by [5]. Local and non local boundary conditions for the numerical solution of wave problems are reviewed in [6]. For the time dependent Maxwell equations in three dimensions an exact nonreflecting boundary condition is derived by [7]. For the numerical solution of linear magnetostatic problems consisting of a homogeneous magnet embedded in free space, the boundary element method (BEM) is ideally suited. Only the surfaces of the magnet need to be discretized and the surrounding space is captured exactly. If nonlinear material behavior is to be considered, the BEM can be coupled with the FEM, see [8]. Alternative (comparable) hybrid methods based on FEM-BEM couplings have been developed for time-harmonic eddy current problems in three-dimensional unbounded domains, see [9]. Standard coupled FEM-BEM often leads to unacceptable numerical costs for a large number of degrees of freedom. Hertel et al have shown in [10] that the memory requirements of these algorithms can be drastically reduced if special hierarchical matrices are used, see also [11]. A comparative study of finite difference based fast Fourier transform methods, tensor-grid methods and the finite-element method with shell transformation for the stray field calculation in micromagnetics is given by [12]. There, the authors investigate the above algorithms with respect to complexity, storage requirements and accuracy obtained for benchmark problems. Depending on the task, one or the other method proves to be advantageous.
Another discretization scheme for the solution of problems in infinite domains is the scaled boundary finite element method (SBFEM), which goes back to [13] and [14], see also [15] and the monograph [16] for a status quo of the SBFEM for bounded domains. This semi-analytical scheme retains advantages of both the finite element method and the boundary element method. In particular, its semi-analytical nature facilitates the accurate representation of radiation damping when modelling acoustic or dynamic problems in unbounded domains, see e.g. [17] and [18], respectively. An application of the SBFEM to two-dimensional electromagnetic field problems can be found in [19]. In a 3D context the surface of the domain of interest is discretized using 2D finite element meshes. It should be noted that the semi-analytical treatment of unbounded regions using SBFEM requires a certain piecewise regularity of the domain.
Our proposed concept is formally placed in the class of truncation of outer boundaries. Here we apply a static condensation, which can be interpreted as a model order reduction technique. We proceed as follows: First, we define the boundary of the outer domain. The definition of the outer boundary of the region under consideration is based on the assumption that here the potential (i.e. the Dirichlet boundary condition) or the normal derivative of the potential (Neumann boundary condition) tends towards zero.
Second, we impose another (inner) boundary separating the solid from the outer space. This would provide the reduced system with the smallest number of degrees of freedom. However, we can also draw an arbitrary boundary within the region between the outer boundary and the boundary containing the solid. The only condition is that the region (denoted as exterior or outer domain in the following) between the inner and outer boundary behaves linearly.
Third, we perform a static condensation of the discretized exterior domain onto the degrees of freedom of the inner boundary.
Finally, a nonlinear finite element analysis of the condensed problem is performed. Since, as mentioned above, the assumption of homogeneous boundary conditions is essential for this procedure, a substitution of the basic variable, here the magnetic potential, is necessary. The solids to be investigated are placed in an externally imposed, spatially constant magnetic H-field. For the external domain, this results in a linear distribution of the potential on the boundary. We now split the magnetic field, which is determined later, into the constant applied field and a fluctuation part. Thus, the descriptive differential equation is formulated with respect to the fluctuations, which by definition disappear on the outer boundary. Furthermore, the challenge arises that the imposed external field is variable, which requires an update of the condensed right-hand side. The numerical implementation of the proposed procedure is done within the finite element framework of AceGen and AceFEM, compare [20,21].

Definition of the boundary value problem
This section presents the magnetic quantities of interest and the equations to be solved. For an introduction to magnetism and magnetic materials, refer to standard textbooks, such as [22,23] or [24]. Let B be the domain of interest, parametrized in x, which is subdivided into an interior B int and an exterior domain B ext . To illustrate the problem, Fig. 1 shows a circular region with a circular inclusion. The outer domain is represented by the light grey circular ring and the inner domain by the yellow and orange inclusion. The interface between the interior and exterior domain is denoted by ∂B b and the boundary of the outer domain is ∂B ext . Let us now consider  The discretization of the boundary of the interior domain ∂B b ist denoted by ∂B int with the outward unit normal n int . In analogy we write ∂B ext,i for the discretization of the same interface associated with B ext with the outward unit normal n ext =: n. At associated points of both discretizations of ∂B b the relation n int = −n ext holds. Finally, the discrete representation of the outer boundary is ∂B ext,a .
As a generic example we consider a boundary value problem of magnetostatics. The Gauss law of magnetostatics is where B is the magnetic induction. Let ϕ denote the magnetic potential. Then, the magnetic field vector is defined as Possible boundary conditions are the specification of a magnetic potential ϕ 0 and a magnetic flux ζ 0 : The weak form of (1) appears as with the virtual magnetic potential δϕ and the virtual magnetic field By definition we assume linear behavior in the exterior domain. In this region we use the constitutive relation where μ 0 denotes the constant magnetic permeability. Now we split the weak form into the two terms where we substitute the constitutive linear relation (7) in a tensorial manner μ = μ 0 1. We start our considerations with an externally applied, constant magnetic H-field over the whole domain B. If there is no magnetic body with a different permeability located within B, no interference fields (stray fields) can occur. If a magnetic body B sol is placed in the external H-field, a magnetization will be generated within the solid and it holds B = μ 0 (H + M). Now, sources and sinks of the magnetization appear on the surface of the magnet, but according to Maxwell's equation they must not exist, since div B = 0. Thus, sources (sinks) of the magnetization must be understood as sinks (sources) of the magnetic field strength. Mathematically, this can be expressed as follows This creates a magnetic field that opposes the magnetization inside the magnet, the demagnetizing field, also known as the stray field outside the magnet. Compare [24,25] for further insights into the topic. In the following analysis we are interested in the computation of magnetization fields (M-field) under consideration of the magnetic stray field (demagnetization field). The stray field in the external environment outside the magnet, is the magnetic field (H-field) generated by the magnetization M in a magnet. The typical task is now to calculate the above fields due to an externally applied constant H-field for a magnet in a large external space. This external space represents, for example, a vacuum or the air environment. The stray fields distort the applied constant Hfield, especially in the immediate vicinity of the magnet, and decay from a certain distance so that they no longer have any effect on the applied field. Based on these preliminary considerations, an H-field H is applied, which we increase via the parameter λ ∈ [0, 1]. This facilitates the substitution of the magnetic potential as Furthermore, this results in the following substitutions for the virtual magnetic potential and magnetic field With these modifications, the weak forms (8) result in For later use, we can now prescribe the magnetic potential on the complete outer boundary ∂B ext,a ; thus ∂B ext,a = ∂B ext,a ϕ and ∂B ext,a B = ∅. If we take the substitution into account, it follows that the fluctuations on the outer rim must disappear by definition, i.e.
Before we turn to the discrete formulation of the boundary value problem in the next section, we perform the formal linearisation of equation set (12). The linearization LinG ϕ = G| n + G ϕ yields the linear increments

Discrete formulation
By means of discretization of the body B with standard finite elements B e we obtain the representation where the index h characterises the discrete areas or elements. The number of elements in the exterior and in the interior domain are denoted by num ext ele and num int ele , respectively; thus, the total number of elements is num ele = num int ele + num ext ele . Due to the substitution of the magnetic potential by the sum of a basic affine component and a fluctuation field, we only have to discretize the fluctuation field and the virtual fluctuation field for a typical finite element B e : where IB I and IB e denote the node-wise and element-wise Bmatrices, respectively, containing the Cartesian derivatives of the shape functions.

The interior domain B int
Substituting Eq. (17) into the weak form and its linear increment, we obtain with the element right-hand side r e,int and the element stiffness matrix k e,int . If we now assemble the algebraic system of the inner domain, we obtain the representations with the stiffness matrix K int and the right-hand side R int associated with the interior domain. In view of the subsequent reduction of the total system, we partition the system of equations with respect to the degrees of freedom d i inside the interior domain B int \ ∂B int and the degrees of freedom d b on the interphase ∂B int between both domains, see Fig. 1: being the incremental global vector of unknowns.

Discretization of the external domain B ext
Substituting the constant and a fluctuation component of the magnetic field into the linear constitutive equation (7) yields For the exterior system B ext the weak formulation is expressed as Since λ and the given H-field (H and therefore B) are constant in an incremental step, the linear increment results in The algebraic system of equations for the external domain B ext appears in the form with the global matrices, after carrying out the classical assembly procedure for B ext , With the aim of further reducing the size of the system of equations, we now partition the degrees of freedom, analogous to the procedure for the inner domain, according to the degrees of freedom d b on the interface between B ext and B int and the degrees of freedom d a within the outer domain with D representing the global solution vector. Since the inner domain is generally characterised by nonlinear behavior, an iterative procedure must be applied to solve the entire system. Here, it is important to note that we also apply a load increase in the (linear) outer region via the parameter λ and that the right-hand side R ext depends on the current fluctuation field, i.e. R ext (λ,D).

Static condensation of the exterior domain
In a next step, we eliminate from equation (26) the degrees of freedom D a , where we have to take into account the explicit expressions for the right-hand side. Starting from we first eliminate D a from (27) 2 : Since the symmetric matrices K ext aa , K ext bb and the matrices K ext ab = (K ext ab ) T are constant, it is advisable from an algorithmic point of view to introduce the auxiliary matrix Substituting Eqs. (28) and (29), we obtain with the constant matrices K ∞ bb and R ∞ b . Finally, we must carry out the assembling procedere of the discrete inner and outer domain, i.e. the assembly of the algebraic systems (30) and (20). This yields the final representation The algorithmic steps are summarized in Table 1.

A simplified nonlinear magnetic model
To verify the performance of the proposed numerical scheme, where the degrees of freedom of the external domain have been eliminated by static condensation, we apply a simplified nonlinear constitutive law for the solid inside B int for the magnetic induction B =B(H) as a function of the magnetic field H. In general, this relation can be expressed via the magnetic permeability μ and thus via the magnetic susceptibility χ : where M denotes the magnetization. The magnetic susceptibility is not constant in ferromagnets, but a non-linear function of the applied magnetic field strength and of the history of the material, i.e. χ =χ (M, ...). In ferromagnetic materials, a so-called saturation magnetization is generally observed, i.e. a material-dependent maximum value for the magnetization is reached. If we consider χ as differential susceptibility, i.e. χ = ∂ H M, then this becomes zero when the saturation magnetization is reached. For the sake of simplicity, we will limit ourselves here to a magnetization curve without hysteresis. This is a characteristic of soft magnets, which do not show much hysteresis. However, the proposed procedure is also valid for the simulation of hard magnets with distinct hysteresis properties. Of course, this model also shows saturation. As a prototype for this we use the function otherwise we obtain Note: This is a very simplified model to describe the magnetization in a solid. Classically, one assumes the existence of a free energy function U (m 1 , m 2 , m 3 ), which is formulated in series expansion of the direction cosines (m i |i = 1, 2, 3) of the magnetization relative to the crystal axes. In our consideration we start from an enthalpy function H(H) = U (m) − μ 0 M S m · H. For our chosen model the  Fig. 2 The normalized response curve for the load case described above enthalpy function for the magnetization part is For the discussion of phenomenological models of hysteresis we refer to [26,27]. To demonstrate the feasibility of the model presented above a homogeneous boundary value problem is considered. The domain of interest corresponds to a square with dimensions of 4 μm × 4 μm and the material parameters as defined in Table 2. In order to generate a magnetic response, the area is subjected to the external magnetic fieldH = λ 0.07 [1, 0] T A/ μm, where λ is a load factor alternating between minus one and one. The corresponding response function can be found in Fig. 2, where · denotes the Euclidean norm.

Stray field calculation with heterogeneous exterior domain
The method presented above is suitable for the approximation of large outer domains, while other established methods, such as the boundary element method or the scaled boundary finite element method, provide solutions for infinite domains. While the established methods can only give solutions for outer domains with a certain regularity, the proposed method has the decisive advantage of being able to take heterogeneous exterior domains into account. This means that interference fields in the outer space, e.g. generated by adjacent or nearby magnetic bodies, can also be taken into account. This is demonstrated in the following example. Three magnetic inclusions, identified as B ext 1 (green), B ext 2 (blue) and B int 2 (red), are considered in direct proximity to one another. In the following, all distances are related to the coordinate system shown in Fig. 3. The green inclusion B ext 1 is described by the center point P 1 (-80 / -110) μm and the radius r 1 = 35 μm. The blue inclusion is also characterised by a centre point P 2 (80 / -130) μm and the radii r 2 = 70 μm and r 3 = 30 μm. The dimensions of the remaining coloured domains can be found in the drawing in Fig. 3. All three inclusions are surrounded by free space, here indicated by B ext 3 (light blue) and B int 1 (orange). To generate a magnetic response the whole domain is treated with an external magnetic field H 2 = 6.3 ·10 −3 A/ μm.
Since the magnetic particles are very closely located, their magnetic fields will influence each other. Due to the fact, that the magnetic stray field flattens out towards zero for large distances to the magnetic particles, a huge airspace must be considered. In this contribution, the airspace is taken into account in two different ways. Firstly, a reference solution was generated using the truncation method. Therefore, the FEM is applied to approximate the whole boundary value problem as sketched in Fig. 3, i.e. the airspace and the magnetic inclusions. Secondly, the static condensation is applied to the outer domain B ext = B ext 1 ∪ B ext 2 ∪ B ext 3 of the boundary value problem given in Fig. 3. This means that the outer domain is reduced onto the border of the interior domain B int = B int 1 ∪ B int 2 , so that no further discretization of B ext is required. Hence, the area to simulate is reduced to B int . The aim is to verify that both results match and to evaluate the performance of both methods. For the following calculations, identical magnetic properties are assigned to the magnetic inclusions, i.e. μ mat = 1.25 · 10 5 (ng μm)/ (A 2 μs 2 ). The free space is considered to be vacuum with μ air = μ 0 being the magnetic vacuum permeability. The material behavior is approximated by a linear material model. Different numbers of degrees of freedom are assigned to the inner and outer domain. The much smaller inner domain has almost 9,000 degrees of freedom, while the outer domain has over 70,000 degrees of freedom. This corresponds to a total number of approximately 80,000 degrees of freedom in the overall system. Of course, the reference simulation provides results for the entire area, but since the outer space is very large, only a fraction of the boundary value problem is used for graphical representation. The cut out used for graphic processing includes the magnetic inhomogeneities and is shown below in Fig. 4. Due to the close position of the inclusions to each other, the field lines corresponding to H show a strong swirling of the fluctuation field around the inclusions, as shown in Fig.  5a. In contrast to the fluctuation field H, the external magnetic field H is a spatially constant distributed quantity. It has the same amplitude in every spatial point of the observed area B = B ext ∪ B int . The external magnetic potential ϕ delivers a linear distribution over B. For this reason, the external magnetic field is not shown graphically. Since the fluctuation field evolves cloud-like and the external field is constant, the entire magnetic field H = H + H also changes spatially. The magnetic field H is derived from the entire potential ϕ = ϕ + ϕ. Hence, Fig. 5b presents the potential ϕ as a coloured contour plot and the derived magnetic field H as streamlines.
The results of the reduced simulation must be identical to those of the reference solution within the internal domain B int . Since the outer space B ext has already been taken into account in the preliminary static condensation, the domain to discretize with finite elements reduces to B int , compare Fig. 6. Therefore, the plots in Fig. 6 c and d only show the evolution of the magnetic potentials ( ϕ and ϕ) and the corresponding fields ( H and H) for the inner domain B int . Similar to Fig. 5, the corresponding potentials are shown as coloured contour plots and the derived magnetic fields are superimposed as stream plots. It is obvious that these variables behave analogously to those from the simulation that was calculated previously. The magnetic inclusion B int 2 coloured in red in Fig. 6 a and b is responsible for most of the turbulence in the illustrations. However, the influences of the adjacent inclusions (B ext 2 and B ext 3 ) are also clearly visible, since they were previously taken into account by the static condensation. The comparison of the potentials and fields presented in Figs. 5 and 6 already confirms a very good agreement of the observed quantities within the reduced domain. In order to emphasize the accuracy of the proposed method, discrete values of the reference solution (whole domain) and the solution obtained on the boundary value problem defined by only the interior domain, are taken along the marked intersection lines A-A and B-B (compare Fig. 6), see Fig. 7. The plots of the discrete values below are congruent and confirm the accuracy of the proposed scheme. A characteristic of classic FEM stiffness matrices is a sparsely populated band structure. This is in complete contrast to the stiffness matrices of the reduced systems, which are typically fully or very densely populated matrices. In Fig.  8 both types of matrices are shown. In order to emphasize the differences between these two matrix types, the sparsely populated matrix of the pure FEM calculation is shown together with the significantly smaller, but dense, reduced matrix in Fig. 8a. The reduced matrix can also be seen in scaled size in Fig. 8b.
In order to evaluate the performance of the presented method against the reference simulation, the times required for solving the system of equations are compared. All following calculations are done using a standard laptop with an i5 processor and 16 GB of RAM. The system of equations described above was reduced in less than ten seconds. Since the reduction of the outer space was carried out as a  (Fig. 6) are given as circular points Fig. 8 The population of the stiffness matrix corresponding to the full FEM solution is presented in a) together with the stiffness matrix corresponding to the reduced calculation marked by a red frame. The population of the reduced stiffness matrix is also presented in b) as a scaled plot preliminary calculation, it does not affect the actual solution procedure of the finite element simulation. Hence, the time needed for the reduction procedure is not included in the timing. Time savings are particularly noticeable in dynamic calculations, therefore the magnetic field applied to the boundary value problem described above is increased from zero to H 2 within 100 load steps. Since the system of equations is solved repeatedly, the time difference adds up. The results obtained for both the reduced and the full FEM simulation are presented in Table 3. A large difference in time and therefore a significant gain of efficiency is obtained.
Although the performance advantage can already be seen for the linear simulations, the reduction will probably be of even greater advantage, especially for nonlinear and timedependent calculations that require the consideration of an external area. In the case of a nonlinear simulation, several Newton iterations are usually required to solve the existing system of equations. Since each individual iteration takes longer for large systems of equations compared to reduced ones, the time savings for serial evaluations of systems of equations, as it is the case with time-dependent problems,  are obvious. At this point, the same boundary value problem as utilized for the linear case is used, but the nonlinear material behavior, which is introduced in Sect. 4, is considered in B int 3 . Similar to the previous simulations, a reference solution is created with the truncation method, which is compared to the reduced method afterwards. For this purpose, the external magnetic field H 2 is increased within 100 load steps in the y direction. The required times of the truncation and the reduced method for solving the systems of equations serially are given in Table 4. has an approximated radius of 70 μm. In this example, the domains B ext and B int 1 are considered to represent vacuum, while a permeability of μ mat = 1.25 · 10 5 (ng μm)/ (A 2 μs 2 ) is assigned to the magnetic inclusion. The exact geometry, including the dimensions, can be found in Fig. 9. The radii of the extended outer domain B ext , as well as of the inner domain B int are presented in Table 5 with the corresponding numbers of degrees of freedom within the domains and on the surfaces. By means of this comparison, the size difference between the two systems is emphasised. Because a single inclusion is considered to be placed within a free space of vacuum, no adjacent sources of interference are obtained. Therefore, the externally applied field H only generates a reaction within the magnetic inclusion B int 2 , so that the magnetic potentials ( ϕ and ϕ) and fields ( H and H) can propagate without being disturbed. Since the results of the static condensation show good agreement with those of the truncation method, only the field curves in the reduced area B int are shown below. The fluctuation potential shown in Fig. 10a shows the typical cloud-like curves. Hence, the fluctuation field H in Fig. 10b results in a spatially nonlinear distribution. Thus, the overall magnetic field H in Fig. 10c also results in a spatially nonlinear distribution, which is particularly pronounced within the area of the inhomogeneity.
To determine the time gain of a dynamic simulation through static condensation, the magnetic field is increased within 100 load steps, from zero to the maximum field strength of H 2 . Both simulations are timed and the results are compared in Table 6. The time advantage of the reduced   . 11 The population of the stiffness matrix corresponding to the full FEM solution is presented in a) together with the stiffness matrix corresponding to the reduced calculation marked by a red frame. The population of the reduced stiffness matrix is also presented in b) as a scaled plot simulation with a factor of 4.66 compared to the full FEM calculation is obvious. For larger systems an increase in time for the static condensation can be seen. However, the time spent for this calculation can be made up quickly, especially with dynamic simulations. That means, a reduction in three dimensions can also be worthwhile, if results in the outer space are not necessarily required.
Analogously to the two-dimensional FEM, the threedimensional FEM also results in sparsely populated stiffness matrices. Since the reduced matrix is very densely populated, this advantage is lost through the static condensation of the outer domain B ext onto the boundary of the interior domain ∂ B int as in the two-dimensional case. However, this disadvantage will be counterbalanced by the advantage of a significantly smaller reduced stiffness matrix. To illustrate the size differences, a comparison of the two matrices can be found in Fig. 11.

Conclusion and outlook
In the present work a method for the consideration of large, heterogeneous, unbounded domains in nonlinear finite element analyses of magnetostatic problems was presented. Compared to established methods such as the BEM or the SBFEM, this way of reducing the external space allows us to take interferences outside the reduced area into account. Due to the static condensation of the outer space, the system stiffness matrix is no longer sparsely populated, but very dense. However, with clever reduction it can be orders of magnitude smaller than the stiffness matrix of the pure FEM. Several numerical examples were used to illustrate the advantages of the statically condensed systems over a full FEM discretization of the exterior domain. For both linear and nonlinear material behavior in the reduced area, large differences in the required computing times of the proposed method and a pure FEM simulation have been observed. Since the method presented is a general scheme for representing external spaces, it can also be applied to other subject areas that are not necessarily related to magnetism.