Efficient Reduced Magnetic Vector Potential Formulation for the Magnetic Field Simulation of Accelerator Magnets

The major advantage of reduced magnetic vector potential formulations (RMVPs) is that complicated coil structures do not need to be resolved by a computational mesh. Instead, they are modeled by thin wires, whose source field is included into the simulation model along Biot-Savart’s law. Such an approach has already been successfully employed in ROXIE for the simulation of superconducting Large Hadron Collider magnets at CERN. This work presents an updated RMVP approach, which significantly outperforms the original method. The updated formulation is postulated, implemented, verified, compared to the original formulation, and applied for the simulation of a quadrupole magnet. The promising results of this work encourage further investigation towards an updated simulation framework for next-generation accelerator magnets.


I. INTRODUCTION
H IGH-TEMPERATURE superconducting (HTS) technol- ogy is expected to have a significant impact on nextgeneration synchrotrons [1].With this step forward, however, magnet design is confronted with new challenges regarding the design of large, high-field and high-quality magnet systems.Computer-aided design and numerical field simulation generally play a crucial role in designing and optimizing accelerator magnet systems, and will do even more so regarding the next-generation HTS magnet systems.For the past decades, the simulation software ROXIE [2] proved to be an indispensable workhorse for designing the low-temperature superconducting (LTS) magnets of the Large Hadron Collider (LHC).ROXIE combines a hybrid finite-element (FE) boundary-element method with a reduced magnetic vector potential (RMVP) formulation [3], [4], leading to fast and accurate simulations [5], [6].Herein, the coils are modeled as thin wires, and their excitation is included as a source magnetic field, which is calculated by Biot-Savart's law.The major advantage of this formulation is that these wires do not need to be resolved by a computational mesh.Especially superconducting accelerator magnets, which typically contain Manuscript received Month Day, Year; revised Month Day, Year.Corresponding author: Laura A. M. D'Angelo (email: dangelo@temf.tudarmstadt.de).
hundreds or thousands of coil windings, greatly benefit from this approach.
Nonetheless, ROXIE and commercial out-of-the-box simulation tools struggle with the multi-scale nature that is imposed by HTS tapes, resulting into excessive computation times [7].The goal of this work is to improve the computational efficiency of the RMVP approach and to contribute towards suitable simulation tools for future HTS magnet design campaigns.To do so, ROXIE's original formulation is first adapted to a pure FE model in Sec.II, and then reformulated in order to diminish the number of Biot-Savart integrals to be calculated to quantify the source magnetic field.This leads to a multistep calculation procedure, which is presented in Sec.III and demonstrated for a case study with an eccentric line current in an iron tube.In Sec.IV, the updated RMVP formulation is analyzed regarding accuracy and performance.Herein, the proposed procedure proves to be clearly superior to ROXIE's original formulation.Sec.V finally showcases the updated RMVP formulation by applying it to the two-dimensional (2D) nonlinear magnetostatic simulation of the LHC's LTS MQXA quadrupole magnet [8].All simulations have been carried out using the freely available FE solver GetDP [9].

II. ORIGINAL RMVP FORMULATION
In this section, the original RMVP formulation from [3] is recapitulated.The original physical problem that has to be solved is the magnetostatic problem with a homogeneous Dirichlet boundary condition, reading Herein, ⃗ A is the sought-for magnetic vector potential (MVP), ν is the reluctivity, which is possibly nonlinear, ν = ν(⃗ r, ⃗ B(⃗ r)), and ⃗ J is the current density, which represents the excitation in this problem.The computational domain V = V a ∪ V i consists of the coil and air domain V a and the iron domain V i , and ∂V is its boundary.Classically, the current excitation is modeled in the right-hand side ⃗ J, e.g. by using winding functions [10].This procedure requires the explicit discretization of the individual wires (or at least half-turns) in the FE mesh.

arXiv:2309.02004v2 [cs.CE] 14 Dec 2023
In contrast, the RMVP method represents the wires by onedimensional (1D) curves, which are not necessarily taken into account in the FE mesh.This benefits the meshing workload in the overall simulation process.The magnetic vector potential (MVP) is decomposed into where ⃗ A s is called the source MVP, and ⃗ A r the reduced MVP.⃗ A s is obtained by evaluating Biot-Savart's law [11] for all spatial coordinates ⃗ r ∈ V .The source domain L ′ represents the line, on which the line current I is located.
In a three-dimensional (3D) model, this would be an arbitrarily complicated 1D curve loop in space; while in a twodimensional (2D) setting, L ′ is reduced to a zero-dimensional (0D) point, which represents a line current going in or out of plane.Furthermore, in 2D, the MVP is assumed to have only a z-component, ⃗ A = A z (x, y)⃗ e z , and Biot-Savart's law (3) becomes Multiple sources are taken into account by superposition.Eventually, the discrete source MVP can be computed on the mesh edges j = 1, . . ., N edge in two ways: One can utilize the partition-of-unity property and calculate the discrete coefficients ⌢ a s,j per edge e j directly by weighting (3) with the j-th edge function ⃗ w j and integrating over that edge e j , ⌢ a s,j = µ 0 4π Alternatively, one performs a weak L 2 -projection of the Biot-Savart integral onto ⃗ A s , This work uses the built-in L 2 -projection of GetDP.The reduced MVP ⃗ A r is computed by solving the boundary value problem (BVP) Using a Ritz-Galerkin approach, the weak formulation is obtained as: Find ⃗ A r ∈ H r (curl; V ), s.t.
is chosen in order to fulfill (9c).Herein, is the tangential trace operator w.r.t. a boundary B [12].The weak formulation (10) is solved by a FE method employing standard edge shape functions.Lastly, the total MVP ⃗ A is composed of ⃗ A s and ⃗ A r following (2).Note that the Biot-Savart integral (3) must be evaluated in the whole domain V , whether one is actually interested in the solution of the whole domain or of only a small sub-domain.

A. Idea and derivation
The domain V is decomposed into a non-permeable subdomain V a (consisting of e.g.air) and a source-free subdomain V i (containing e.g.iron).Thereby, V a is fully enclosed by V i , and Γ denotes the interface between those domains.This configuration is illustrated for a simply 2D case in Fig. 1.The BVP ( 9) is expressed for both domains separately and, after introduction of the source MVP, reads where ⃗ A a and ⃗ A i are the MVPs and ⃗ H a and ⃗ H i the magnetic field strengths in V a and V i , respectively.Then, an additional field ⃗ A m is added to the total MVP in V a , which is the solution of a so-called image problem as will be described in Sec.III-C, such that homogeneous boundary conditions are enforced at Γ, yielding the BVP Hence, the sub-domain solutions ⃗ A b and ⃗ A i are tangentially continuous at Γ. Their co-normal derivatives (magnetic field strengths), however, feature a jump prescribed by ⃗ H s + ⃗ H m .This jump can be interpreted as a surface current density in A/m, By substituting 14) can be transferred to a single-domain formulation w.r.t. the so-called reaction MVP ⃗ A g , in which ⃗ K g is imposed on the interface Γ.This so-called reaction subproblem will be introduced in detail in Sec.III-D.This approach resembles the sub-domain FEM based on a perturbation technique, where a missing continuity at the interface of two sub-domains is also restored via a jump in the tangential magnetic field strength expressed as a surface current density source [13].

B. Ansatz and Biot-Savart sub-formulation
As a result of the previous section, the total MVP ⃗ A is decomposed into where ⃗ A s is called the source MVP, ⃗ A m the image MVP, and ⃗ A g the reaction MVP.⃗ A s is obtained by evaluating Biot-Savart's law (3) for all spatial coordinates ⃗ r ∈ V a that are of interest, but at least for ⃗ r ∈ Γ.This is the major advantage to the original RMVP approach [3], where ⃗ A s needed to be calculated in the whole domain V .

C. Image sub-formulation
The image MVP ⃗ A m is the solution of the BVP This sub-formulation originates from using the method of images [11].In this context, ⃗ A m represents the field that needs to be added to ⃗ A s in V a in order to enforce electric boundary conditions at Γ.In other words, the MVP subtotal ⃗ A s + ⃗ A m is the equivalent of a Green's function obeying a homogeneous Dirichlet boundary condition at Γ [11].
Using a Ritz-Galerkin approach and ⃗ n × ⃗ H m for the tangential component of the image magnetic field strength, the weak formulation is obtained as: A ′ m and ⃗ n× ⃗ H ′ m are corresponding test functions, and H −1/2 (curl; Γ) is a trace space [12].Here, the boundary condition (17b) is weakly imposed by (18b), yielding a saddle-point problem [14].In this way, the quantity ⃗ n × ⃗ H m , which will be needed for the calculation of ⃗ A g , is already at hand.The weak sub-formulation (18) is eventually solved by a FE method employing standard edge shape functions.

D. Reaction sub-formulation
The reaction MVP ⃗ A g is obtained by solving the BVP with ⃗ J g = ⃗ K g δ Γ , where ⃗ K g is the surface current density (15) and δ Γ the delta distribution function defined by Here, ⃗ are the source and image magnetic field strengths, respectively.Visually, the current excitation in V a has been shifted onto the interface surface Γ.The weak formulation of ( 19) reads: is chosen in order to fulfill (19b).The weak subformulation ( 21) is solved by a FE method employing standard edge shape functions.Finally, the total MVP ⃗ A is composed of ⃗ A g , ⃗ A m and ⃗ A s following ( 16).

E. Treatment of nonlinearities
Note that ⃗ A g corresponds to the full MVP in the domain V i .This facilitates the treatment of nonlinearities, as only ⃗ A g is affected by a nonlinear reluctivity ν = ν(⃗ r, ⃗ B(⃗ r)).To approximate the nonlinear material characteristic, a standard nonlinear iteration scheme can be employed such as fix-point iteration or Newton's method [15], which is then only applied to the reaction sub-problem (21).Other than that, no further modifications are needed to take nonlinearities in the updated RMVP formulation into account.

F. Illustration: Eccentric line current in an iron tube
Consider as a case study an infinitely long eccentric line current I in air surrounded by an infinitely long iron tube.This model's 2D cross-section is shown in Fig. 1.V a is defined as the air and coil region within the iron tube, while the iron tube itself is chosen as V i .Thus, the interface Γ is the circular boundary between those two domains.
The RMVP formulation is applied to compute the MVP and magnetic flux density caused by the direct current I. Fig. 2a shows the source MVP ⃗ A s , which is obtained by evaluating the Biot-Savart integral (3).It represents the MVP as if the eccentric line current were located in free space.The corresponding image MVP ⃗ A m calculated by ( 18) is

A. Benchmark model: 2D racetrack coil
The RMVP formulation is implemented in the freely available open-source FE solver GetDP [9] and employed to carry out a 2D linear magnetostatic simulation of a racetrack coil in an iron yoke.To justify the 2D approximation, it is assumed that the racetrack coil is very long compared to its crosssection diameter, such that the effects at the coil end windings V eval (a) Two coil windings (hatched rectangles, see Fig. 3b for a detailed sketch) in air (white, denoted as Va) surrounded by an iron yoke (gray, denoted as Vi).The interface boundary between iron and air is denoted as Γ.The dashed circular domain represents the evaluation domain for the magnetic energy.
(b) Detailed view of the coil winding configuration for Nx = Ny = 3.In the volumetric case (1), each winding is modeled as a rectangle with a given volumetric current density.In the RMVP case ( 21), each half-turn is represented by a set of points (purple and gray dots) representing line currents in or out of plane, depending on the current orientation.In the simplest case, each half-turn is modeled by a single line current located in the cross-section center (purple dots).can be neglected at a central cross-section of the model.Fig. 3a shows the geometry, which consists of two winding groups (hatched rectangles) containing wires with rectangular crosssection and embedded in an air domain V a (white), which is surrounded by an iron yoke V i (gray).Γ is the interface between V a and V i .For the iron yoke, a linear permeability of µ i = 4000µ 0 has been chosen.Since the line currents come together with singularities, the magnetic energy strives for infinity.Therefore, the magnetic energy considered below is evaluated in a sub-domain, which is the dashed circular domain V eval within V a .
The configuration of each winding group is shown in Fig. 3b.In the 2D reference model, each half-turn of the coil is modeled as a surface (black rectangles in Fig. 3b) with a given surface current density.In the 2D RMVP setting, each half-turn is discretized by a set of points (see Fig. 3b, gray and purple dots), which represent line currents in or out of plane, depending on the current orientation.In the simplest case, a half-turn is modeled by a single line current located in the cross-section center as seen in Fig. 3b (purple dots).
The results of the RMVP approach are verified against reference values obtained by a conventional 2D FE simulation modeling the half-turns as surfaces and solving (1) while using a very fine mesh.The resulting magnetic flux density magnitude of the proposed method is shown in Fig. 4, where the maximal value has been capped to that of the reference solution.Particularly high fields are obtained at the four inner corners of the iron yoke due to the sharp geometry as well as the linear material properties and, of course, at the line currents due to their singular nature.

Mesh length h
Relative error Fig. 5: Linear convergence of the L 2 -error in the evaluation domain (solid, purple) w.r.t. the mesh length for the RMVP formulation.The L 2 -error in the whole domain (dashed, orange) does not converge due to the singularities.

B. Convergence Analysis
For the convergence analysis, the L 2 -error w.r.t. the volumetric reference solution is considered.The line currents have to be excluded from the considered domain, as otherwise the L 2 -error would not converge due to the singularities.Therefore, a sub-domain V eval ⊂ V is chosen, such that L ′ ⊈ V eval as seen in Fig. 3a (circular region).
For lowest order FEs, one would normally expect a quadratic convergence of the L 2 -error.However, this only holds if, among other criteria, the right-hand side of the BVP is in L 2 (V ) [16].In the reaction sub-problem (19), the righthand side is a Dirac source term, making the problem not regular.Although classical convergence results are therefore invalid here, it has been shown that the L 2 -error of 2D elliptic problems with Dirac source terms converges linearly [17].Indeed, this is observed for the L 2 (V eval )-error in Fig. 5

(solid purple line).
Varying the size of V eval , and thereby the remoteness of ∂V eval to the line currents, neither improves nor worsens the linear convergence behavior, as long as the V eval excludes the line currents or their immediate neighborhood.Then, the magnetic energy gets highly overestimated and does not converge at all (see Fig. 5, dashed orange line).

C. Performance comparison to the original formulation
The runtimes of the original and updated RMVP formulations are compared for the 2D magnetostatic linear simulation of the benchmark model with more than 100, 000 degrees of freedom on a standard workstation.Both the total runtimes and the runtimes of only the Biot-Savart integral evaluation are measured.Additionally, the updated RMVP simulation is carried out for the optimal case that only the source MVP on Γ has to be computed as well as the worst case, in which the source MVP in V a is wanted.The latter would be the case if one is interested in the magnetic field in the whole aperture of a magnet.In practice, magnet designers are often interested in field values at particular points, e.g. in a circular curve for a multipole coefficient analysis to investigate the field quality of the magnet [18].Therefore, the average runtime is expected to range between those two extreme cases.
Figure 6 shows the measured total runtimes (purple) and the portions of the Biot-Savart law evaluation therein (orange).The following two observations are made, from which two important conclusions are drawn: 1) The Biot-Savart integral computation heavily dominates the total runtime of both RMVP formulations.Therefore, improving that computation will significantly enhance the whole RMVP procedure.In favor of the RMVP formulation, the source field computation has been a research topic for decades and several efficient solution techniques exist to improve runtime, e.g., by utilizing a reduced scalar potential [19], [20], exploiting the existence of closed form expressions for conductors of specific geometric shapes [21], or employing fastmultipole methods [22].Last but not least, one could also parallelize the computation of the Biot-Savart integrals.
2) The updated RMVP formulation is by far computationally superior to the original formulation.This holds true even for the worst case.It is even more obvious when used for the calculation of the magnetic field at particular points, along certain curves or in small-sized regions in the magnet's aperture.Then, the updated method impressively outshines the original procedure.Hence, Fig. 6 illustrates the performance gain of the updated RMVP formulation with respect to the standard formulation and at the same time, parallelization of the Biot-Savart solver as a straightforward and promising measure for further improvement.

D. Distance of the line currents to the interface boundary
The numerical quality of the updated RMVP formulation significantly depends on the distance of the line currents to Γ, which typically coincides with the interface to the iron yoke, but does not have to.To study that behavior, an eccentric line current in air with homogeneous Dirichlet boundary conditions is considered.For this problem, an analytical solution exists and will serve as a reference.While the position of the line current is fixed, the radius of Γ, R Γ , which is an artificial rather than a material interface, is changed.Let δ denote the distance of the line current to Γ.The L 2 -error of the However, the error rapidly decreases with increasing distance δ: Already a tiny gap of δ ≥ 10 −3 m or ∆ → 10 −1 between Γ and the line current is sufficient to eliminate this source of error completely.Since it is realistic for coils and current sources in general to be positioned a few millimeters away from the iron yoke, this numerical behavior is expected to be negligible in practical applications or at least easy to mitigate by substituting the outer coils by line currents which keep sufficient distance to Γ.

E. Application to three dimensions: 3D racetrack coil
The RMVP formulation as presented in Sec.III can be applied to a 3D model with only slight modifications compared to the 2D procedure: 1) The interface Γ between V a and V i is now represented by a 2D surface instead of a 1D curve.
2) The coils are now being modeled by 1D current lines instead of 0D points.Accordingly, the general Biot-Savart expression (3) has to be utilized.For the numerical integration, the current lines are approximated using discrete line elements.
3) The MVP can now depend on all three space coordinates and can point in any direction, instead of depending  only on x and y and featuring only a z-component as in the 2D approximation.Appropriately, 3D edge shape functions have to be used for the FE discretization of the 3D MVP. 4) The MVP is not intrinsically gauged as it is the case for 2D problems [18].Therefore, an explicit gauging has to be introduced in order to enforce a unique solution [23].
Here, a Coulomb gauge [18] is employed, but one could also use a different gauging technique such as a treecotree gauge [24].Apart from those aspects, the core RMVP formulation does not change for 3D problems.
The 3D RMVP formulation is demonstrated on a 3D racetrack coil model with an elliptical iron yoke and nine windings represented by nine 1D closed curves as shown in Fig. 8.As for the 2D racetrack coil problem in Sec.IV-A, a linear iron permeability of µ i = 4000µ 0 is considered.The resulting magnetic flux density is shown in Fig. 9 and Fig. 10 in the whole model and on the central xy-cross-section, respectively.As in the 2D case, particularly high fields (recognizable as yellow arrows) are obtained near the line currents, where the field strives to infinity.
FiQuS is an open-source, free Python package, which enables a sustainable, consistent and reproducible workflow for computational engineers developing accelerator magnet models.This is especially important in the area of particle accelerator design, where the development of accelerator components spans over generations of researchers and engineers of multiple disciplines [27].The updated RMVP formulation is implemented into FiQuS as a generalized template, which FiQuS adapts during runtime according to user-defined geometrical and physical information about the model.
Fig. 11 shows one-eighth of the geometry of the considered MQXA quadrupole, which contains 61 windings (white rectangles), resulting into 488 windings in total [25].Trapezoidal copper spacers (hatched trapezoids) help to achieve the desired coil form.An iron yoke (plaid area) surrounds the coil configuration.
Each winding is approximated by one line current in its center, leading to 488 line currents and individual Biot-Savart integrals in the source sub-problem in total.A predefined BH-curve from the FiQuS material template database is used to model the iron yoke's nonlinear behavior.The nonlinear  problem is solved by Newton's method, for which GetDP's built-in Jacobian functions are utilized.

B. Simulation results
Fig. 12 shows the computed magnetic flux density in the magnet.The results are compared to a reference simulation performed by a conventional 2D FE method taking the windings into account by surface current densities.Fig. 13 shows the relative difference ϵ r of the magnetic flux density magnitude obtained by the RMVP to that of the reference simulation.As expected, high discrepancies (ϵ r ≥ 1) occur in the neighborhood of the line currents, which represent singularities (clearly visible in Fig. 13 as white dots).Outside the immediate neighborhood of the singularities, a good approximation is achieved with relative differences in the order of ϵ r = 10 −3 , . . ., 10 −1 .An even better alignment between the RMVP and reference simulation is expected if the longish winding shapes were resolved by multiple line currents instead of just one at their centers.
These results verify the procedure of treating nonlinear material characteristics as described in Sec.III-E.Furthermore, they demonstrate the operability of the updated RMVP implementation within FiQuS, with which the formulation could become accessible not only for CERN magnet designers, but also for the global magnet engineering community.

VI. CONCLUSION
This work proposed an updated RMVP formulation for accurate and fast magnetic field simulations of superconducting accelerator magnets.The formulation was postulated and verified against a 2D benchmark model.A runtime comparison showed that the proposed method clearly outperforms the original formulation.However, because of the Dirac source term occurring in one of the sub-problems, the L 2 -error only shows a linear convergence for lowest order FEs instead of a quadratic one.This issue is well understood in the scientific computing community, and different techniques exist to improve the convergence order [28].Furthermore, the updated RMVP formulation was also employed to successfully simulate a 3D racetrack coil model.
Finally, the updated RMVP procedure was embedded in the open-source and free Python package FiQuS, guaranteeing the greatest possible applicability in the magnet engineer community, and successfully employed to carry out a 2D nonlinear magnetostatic simulation of a quadrupole magnet.
The promising results of this work encourage an expansion of that method towards an updated simulation framework for HTS accelerator magnets, including physical formulations suitable for HTS magnets [29] and magnetization models for HTS tapes [30].

Fig. 1 :
Fig. 1: 2D case study model: An infinitely long eccentric line current I in an air domain Va surrounded by an infinitely long iron tube V i .
(a) Equipotential lines of the source MVP ⃗ As in the eccentric wire case study calculated by(3).The source current is depicted as a green line.For the sake of visualization, ⃗ As is calculated in the whole domain instead of only at the interface boundary.(b)Equipotential lines of the image MVP ⃗ Am in the eccentric wire case study calculated by(18).The computation is only done in the sub-domain Va.(c)The green lines visualize the surface current density ⃗ Jg at Γ determined by(15).The equipotential lines of the reaction MVP ⃗ Ag in the eccentric wire case study calculated by(21).(d)Equipotential lines of the total MVP ⃗ A = ⃗ As + ⃗ Am + ⃗ Ag in the eccentric wire case study.The source current is depicted as a green line.

Fig. 2 :
Fig. 2: The MVP components obtained in the sub-formulations of the RMVP formulation and the resulting total MVP for the case study of the eccentric wire in an iron ring.

Fig. 3 :
Fig. 3: Racetrack coil model used for the purpose of numerical studies.

Fig. 4 :
Fig. 4: Magnetic flux density magnitude in the racetrack coil obtained by the updated RMVP formulation.

Fig. 6 : 4 O(∆ − 4 ) 2 - error Fig. 7 :
Fig. 6: Runtime comparison between the original and updated RMVP.Both the total times (purple) and the proportions of the Biot-Savart integral evaluations alone (orange) are depicted.For the updated RMVP case, the worst and best case with computation of ⃗ As in Va and Γ has been considered, respectively.

Fig. 8 :
Fig. 8: Geometry of the 3D racetrack coil model consisting of nine windings, which are represented by 1D closed curves and which are surrounded by an elliptically shaped iron yoke.

Fig. 9 :
Fig. 9: Magnetic flux density in the 3D racetrack coil model, computed by the updated RMVP formulation.The nine windings are modeled by 1D curve loops (black).

Fig. 10 :
Fig. 10: Magnetic flux density at the central xy-cross-section of the 3D racetrack coil, computed by the updated RMVP formulation.The nine windings are indicated in black.

Fig. 12 :
Fig. 12: Magnetic flux density magnitude in the MQXA quadrupole computed by the RMVP formulation.

Fig. 13 :
Fig. 13: Relative difference between the magnetic flux density magnitude computed by the RMVP formulation and of the reference simulation performed by a conventional 2D FE solver.