Numerical simulations for magnetostatics

As Project Elara’s free-electron maser design primarily uses permanent magnets, it is essentially to be able to describe them mathematically. We describe the general concepts in modelling permanent magnets in Guide to modelling permanent magnets. Unfortunately, in many cases, it is not possible to find analytical solutions, and we must use numerical methods instead. This page aims to serve as a basic reference for performing numerical simulations of permanent magnets.

Numerical methods

The typical numerical method used for magnetostatics is finite-element analysis,. To begin, we start with the basic differential equation of permanent magnets¹:

$$\nabla \cdot (\mu (-\nabla {\psi }_m)+\mathbf{M})=0$$

Assuming a constant permeability $\mu$, one therefore obtains:

$$-\mu {\nabla }^2{\psi }_m+\nabla \cdot \mathbf{M}=0$$

To derive the weak form, we multiply by test function $v$, and integrate both sides (see a more in-depth description of the finite element method in RF cavity simulation with aperture). This gives us:

$$\begin{array}{r}{\int }_{\Omega }v(-\mu {\nabla }^2{\psi }_m+\nabla \cdot \mathbf{M})dV=0\\ -\mu {\int }_{\Omega }v{\nabla }^2{\psi }_{m}dV+{\int }_{\Omega }v\nabla \cdot \mathbf{M}=0\end{array}$$

We can simplify this by using integration by parts or by simply consulting a list of vector calculus integral identities. This gives us:

$$\mu {\int }_{\Omega }\nabla v\cdot \nabla {\psi }_{m}dV-\mu {\int }_{\partial \Omega }v\nabla {\psi }_m\cdot d\vec{S}-{\int }_{\Omega }\mathbf{M}\cdot \nabla vdV+{\int }_{\partial \Omega }v\mathbf{M}\cdot d\vec{S}=0$$

Combining the volume and surface terms gives us:

$$\mu {\int }_{\Omega }(\nabla v\cdot \nabla {\psi }_m-\mathbf{M}\cdot \nabla v)dV-\mu {\int }_{\partial \Omega }v(\nabla {\psi }_m-\mathbf{M})\cdot d\vec{S}=0$$

Note that $\mathbf{M}$ must be pre-provided. In our case for permanent magnets (as detailed previously), this is a piecewise function:

$$\mathbf{M}=\{\begin{array}{ll}\frac{B_r}{{\mu }_0}\hat{\mathbf{z}},& \text{inside magnet}\\ 0,& \text{outside magnet}\end{array}$$

Where $B_r$ is the remanence, as before. This can then be solved for suitable boundary conditions (as described above). There are dedicated open-source magnetostatics codes such as Radia as well as general-purpose FEM packages (such as FreeFEM and MFEM) that can be used to solve the resulting BVP numerically.

Numerical magnetostatics in Python

While the finite-element method is the most general and powerful method for magnetostatics simulations, the Python library Magpylib is predominantly used for our calculations of magnetic fields. It implements the analytical formulas for a variety of magnetic field configurations, which have been found for a surprisingly-many number of problems (albeit frequently in terms of special functions or with very complicated solutions)². Therefore, it is capable of modelling permanent magnets, electromagnets, and multipoles (including dipoles) to great accuracy. The linearity of magnetic fields greatly improves the ability to find new solutions from existing ones, and hence the generality of such packages. We have used it to model thin-solenoid magnetic lenses, and will in the future use it to calculate magnetic fields of quadrupole doublet lenses and undulator magnets, which both can be modelled as per their documentation on systems with inhomogeneous magnetization.

Footnotes

1. Referenced from the Wolfram guide on magnetostatic modelling for permanent magnets. ↩

2. See Ortner & Bandeira (2020) ↩