Online magnet calculators

Last week I received the type of tech-help inquiry that accounts for 83% of strokes among software developers. A user reported that results from our 3D code did not agree with numbers he obtained from one of those free magnet calculators that you find on the Internet. Given the choice of believing a finite-element, multi-processor program costing several thousand dollars or a free Web app, which would you choose? Why the app of course! After all, it’s on the Internet!

I wanted to put the issue to rest and to check out some of the calculators, so I set out to resolve the discrepancy. The first step was to add a filter. The user had expended several hours on the problem and, like many tech help requesters, wanted to share the agony. According, he sent a large collection of spreadsheets, Powerpoint files and Magnum scripts. I negotiated things down to the following test problem:

» A magnet in free space has a square cross section (3.5 mm x 3.5 mm), a height of 1.3 mm and a remanence field Br = 0.92 T. The magnetization direction is along the height direction. What is the flux density B at a position 2.5 mm above the magnet (i.e., along the magnetization direction).

The user reported a value 47 mT with an on-line calculator and 37 mT with Magnum.

To begin, here’s a link to an app from IBSMagnet that handles the on-axis fields of circular or rectangular magnets:

http://www.ibsmagnet.com/en/fluxdensity/flussdichte_berechnung.php

Figure 1 shows the setup and the returned value of 47 mT at the test point.

IBSMagnet calculator

Figure 1. IBSMagnet calculator.

As a crosscheck, I found a second calculator from MS-Schramberg at

http://www.magnete.de/en/magnetic-field-calculation/flux-density.html

For this one, you enter one of their materials and it supplies Br. The choice SmCo5 143/143 had Br = 0.92 T. Figure 2 shows the returned value as 44.29 mT.

MS-Schramberg calculator

Figure 2. MS-Schramberg calculator.

I set up my own Magnum calculation rather than employ the user’s setup. For the simple geometry, it took only a few minutes to generate a solution. Here’s a link to a zip archive containing the input files CalcTest.min, CalcTest.gin and CalcTest.scr: CalcTest.zip. Defining the mesh is the critical step. Here is the GLOBAL section from CalcTest.MIN:

GLOBAL
  XMesh
   -10.0  -3.0  0.50
    -3.0   3.0  0.25
     3.0  10.0  0.50
  End
  YMesh
   -10.0  -3.0  0.50
    -3.0   3.0  0.25
     3.0  10.0  0.50
  End
  ZMesh
   -10.0  -3.0  0.50
    -3.0 -0.65  0.20
    -0.65 0.65  0.10
     0.65 3.0   0.20
     3.0  10.0  0.50
  End
  RegName(  1): Air
  RegName(  2): Magnet
END

Besides good resolution for the problem, the mesh has two noteworthy properties:

  • The solution volume is much larger than the magnet and the measurement point to approximate the free-space condition. Variable resolution is used to speed the calculation.
  • The element divisions of the foundation mesh are set up to correspond exactly with the faces of the magnet. In this way, the magnet has precisely the desired volume, even without surface fitting.

The analysis script CalcTest.SCR calls for a calculation at a point on the symmetry axis 2.5 mm from the magnet surface, [0.00 mm, 0.00 mm, 3.15 mm]. The Magnum value is 45.07443 mT, consistent with the online calculators.

There’s no documentation of the methods used by the calculators, but here’s my guess. For the purpose of field calculations, modern magnetic materials like SmCo and NdFe have a wondrous property: the domains are firmly locked along the magnetization direction and are unaffected by other magnets and coils in the vicinity. Therefore, surfaces normal to the magnetization direction carry a uniform surface current density J = Br/μ0. With a known distribution of current, the on-axis flux density can be determined from a simple Biot-Savart integral. For example, a circular disk magnet looks like a finite-length solenoid. The calculation involves a summation over differential current density elements and the accuracy depends on the number of terms. Hence, the small difference between the two calculators.

Finally, I want to emphasize why a code like Magnum costs a lot more than a free calculator. Magnum calculates flux density at all locations (not just the axis), handles assemblies of multiple permanent magnets of any shape and can include the effect of iron (saturated or unsaturated), all with accuracy approaching 1 part in 10^6.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.

Using two-dimensional field solutions in OmniTrak

An OmniTrak feature that may not be familiar is the capability to include two-dimensional field solutions from Estat or PerMag. I recently had a consulting job where this was just the thing I needed, so I thought I would take this opportunity to review why the feature is valuable and how to apply it.

There are two application areas that come to mine where two-dimensional field solutions are useful. The first is when a three-dimensional beam propagates through a symmetric magnetic focusing device (like a solenoidal lens). If we want to include space-charge effects, the electric-field solution is certainly three-dimensional. On the other hand, the beam may have little or no effect on the strong magnetic field. In this case, it is much quicker and more accurate to construct a solution with PerMag rather than Magnum.

The second area can be illustrated with an application example. Figure 1 shows the cross section of a configuration proposed by Harald Enge[1] for an achromatic bending magnet for electrons and ions. The vertical field in the plane z = 0.0 varies as

Bz(x,0,0) = A x^n

Particles enter near the point (0.0, 0.0, 0.0) though a small slot in the left-hand pole at an angle α with respect to the x-axis. The configuration has the property that particles follow a loop trajectory such that they leave the magnet at the same position with angle -α, independent of their kinetic energy. A value n = 0.8 corresponds to a 90° reflection (α = 45°). A value n = 1.0 gives a 81.4° deflection (α = 40.7°). A field with n = 1 is easily achieved — the configuration of Fig. 1 represents one-half of a magnetic quadrupole lens. In this case, the pole faces adjacent to the vacuum gap have a hyperbolic shape, xz = K [2].

PerMag solution for an Enge magnet with n = 1.0

Figure 1. PerMag solution for an Enge magnet with n = 1.0. The red arrows show the magnetization directions in permanent-magnet regions.

 

The bending-field properties hold if the assembly of Fig. 1 has infinite length in y (out of the page). On the other hand, a practical system has finite length. One of the goals of calculations is to determine the effect fringing fields. In other words, adding length in y costs money — what is the shortest length consistent with accuracy? We would like to compare trajectories in a finite quadrupole to those in an infinite system. Years of experience have taught me that it is relatively easy to model real-world devices with three-dimensional codes, but extremely difficult and wasteful to approximate ideal devices like the infinite quadrupole. It’s usually much more effective to use a two-dimensional solution.

For this example, I generated the PerMag field solution of Fig. 1 for input to OmniTrak. The goal was to calculate three-dimensional electron trajectories with horizontal and vertical input displacements in the single-particle limit (TRACK mode). There are two points to note:

  • OmniTrak requires at least one three dimensional mesh to record trajectory vectors. Therefore, it is necessary to include a dummy electric field solution with no applied field. I defined a region covering the volume of anticipated orbits (0.0 ≤ x ≤ 280.0, -100.0 ≤ y ≤ 100.0, -20.0 ≤ z ≤ 20.0) with grounded plates at the top and bottom boundaries in z (i.e., zero electric field everywhere).
  • In planar calculations, the PerMag convention is that the system has variations in x-y with infinite length in z. The field variations in the Enge paper have the same x axis but infinite length in y. Therefore we must rotate the two-dimensional field about the x axis on entry to OmniTrak.

Here is the complete OmniTrak input script to represent electron beams at 2.0, 4.0, 8.0, 16.0 and 20.0 MeV entering at α= 40.7° with a horizontal with of +-3.2 mm:

FIELDS
  EField3D: DummyField.HOU
  BField2D: PMSolution.POU
  Rotate B 90.0 0.0 0.0
  DUnit: 1.0000E+03
END
PARTICLES TRACK
  PLIST
* Mass Charge Ke      x     y    z     ux      uy      ux
* =========================================================
   0.0 -1.0  2.0E6  0.001 -5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  2.0E6  0.001  0.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  2.0E6  0.001  5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  4.0E6  0.001 -5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  4.0E6  0.001  0.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  4.0E6  0.001  5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  8.0E6  0.001 -5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  8.0E6  0.001  0.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0  8.0E6  0.001  5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0 16.0E6  0.001 -5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0 16.0E6  0.001  0.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0 16.0E6  0.001  5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0 20.0E6  0.001 -5.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0 20.0E6  0.001  0.0  0.0  0.7501  0.6613  0.0000
   0.0 -1.0 20.0E6  0.001  5.0  0.0  0.7501  0.6613  0.0000
  END
  DT  1.0E-12
END
DIAGNOSTICS
  PARTLIST
END
ENDFILE

Here, DummyField.HOU is the dummy electric field solution to define the three-dimensional mesh and PMSolution.POU is the field calculated by PerMag. The time step of 1.0^{-12} s corresponds to a vector length of 0.3 mm. For this choice, the particle list shows that energy is conserved to better than 0.5 parts in 10^6. Figure 2 shows the resulting trajectory projections in the plane z = 0.0. The results are not quite perfect because the calculated field for the geometry of Fig. 1 had n = 0.99.

Electron trajectories

Figure 2. Electron trajectories at 2.0, 4.0, 8.0, 16.0 and 20.0 MeV.

Footnotes

[1] H.A. Enge, Achromatic magnetic mirror for ion beams, Rev. Sci. Instrum 35 (1963), 385.

[2] My book Principles of Charged Particle Acceleration discusses magnetic quadrupoles on page 96.

[3] Contact us : techinfo@fieldp.com.

[4] Field Precision home page: www.fieldp.com.

Mapper: creating field maps from midplane data

Information on 3D electric and/or magnetic fields can be ported to the OmniTrak program for charged-particle beam physics in two forms:

  • Finite-element field solutions from HiPhi or Magnum.
  • Maps of field values on a regular grid.

With regard to the second option, the Mapper program supplied with OmniTrak is one source of field maps. In Mapper, the user enters the geometry of poles or electrodes to derive a midplane approximation for dipole-type fields. Approximate off-midplane values are then determined from field expansions.

In a recent consulting job to characterize a spectrometer, I was supplied with a measured set of midplane field values for a dipole magnet. I realized that it would take relatively small modifications to include such data in Mapper calculations. Given midplane values, I could take advantage of existing Mapper features to calculation off-midplane fields and to create maps that could be imported directly into OmniTrak.

The measured field values were contained in a spreadsheet of Bz values on a regular grid of the form

 x0   x1   x2   ...
 y0   B00  B10  B20
 y1   B10  B11  B21
 y2   B20  B21  B23
 ...

This ordering was well-matched to that of OmniTrak field maps. I copied the block of Bij values in Open Office Calc and pasted it into the TextPad text editor (in tab-delimited form). Using the capability to operate on regular expressions, I then did a global search and replacement,  changing \t (tab) to \n (return). The resulting document had entries in the sequential form

 B00
 B10
 B20
 ...
 B10
 B11
 B21
 ...

The values could be read with the loop

DO J=0, JMax
  DO I=0, IMax
   ...
  END DO
END DO

In the standard midplane map file, the data values are preceeded by a header of the form:

 B          (Map type, E or B)
 58         IMax
 15.0       XMin
 305.0      XMax
 40         JMax
 80.0       YMin
 280.0      YMax
 20         KMax
 -10.0      ZMin
 10.0       ZMax
 1000.0     DUnit (Units per meter)
 10000.0    BConv (Units per tesla)

The x-y values in the header define the characteristics of the midplane data. The z values give the desired vertical range of the three-dimensional map. In the example, the value DUnit = 1000.0 implies that coordinates are specified in millimeters and BConv = 10000.0 implies that the listed field values are in gauss.

I added the new command Create table from data to Mapper (Fig. 1). The command initiates a Load file dialog, showing input files with names of the form FileName.DAT. The program handles everything else automatically, creating the output field map FileName.MTX in text format. This file may be ported to OmniTrak, or loaded into Mapper for plotting and editing (Fig. 1).

Mapper screenshot

Figure 1. Mapper screenshot showing a three dimensional field map constructed from midplane measurements with the new Create table from data command.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.

Custom trajectory color-coding in Trak

A customer recently asked about color-coding particle trajectories in Trak and OmniView. He wanted to start particles from an input list (PRT) and to plot trajectories with colors determined by the injection kinetic energy. This is a fairly tall order, but I came up with a procedure for plotting up to three classes of particles using the existing features of the codes.

The strategy is to use the plot filter capability built into Trak and OmniView. Figure 1 shows the filter dialog. Depending on minimum and maximum settings of quantities, model particles may be divided into three classes with trajectory plots in blue, red and green. The problem is that the allowed filtering variables are the model particle current, charge and mass.

Filter dialog in Trak

Figure 1. Filter dialog in Trak.

A simple solution is to introduce small changes of the filtered allowed quantities in the PRT file or Trak script that correlate with the desired properties. To illustrate, here is a sample PLIST input section from a Trak script:

PLIST
*  Mass  Chrg   KEng      x    y      z     px    py    pz
* =========================================================
   0.0  -0.99999  10.0E3  0.00  0.00  -2.499  0.00  0.00  1.00
   0.0  -0.99999  10.0E3  0.05  0.00  -2.499  0.00  0.00  1.00
   ...
   0.0  -0.99999  10.0E3  0.65  0.00  -2.499  0.00  0.00  1.00
   0.0  -0.99999  10.0E3  0.70  0.00  -2.499  0.00  0.00  1.00
* =========================================================
   0.0  -1.00001  12.0E3  0.00  0.00  -2.499  0.00  0.00  1.00
   0.0  -1.00001  12.0E3  0.05  0.00  -2.499  0.00  0.00  1.00
   ...
   0.0  -1.00001  12.0E3  0.65  0.00  -2.499  0.00  0.00  1.00
   0.0  -1.00001  12.0E3  0.70  0.00  -2.499  0.00  0.00  1.00
ENDLIST

Electrons with kinetic energy 10.0 keV have charge q = -0.99999*qe and those with kinetic energy 12.0 keV have q = -1.00001*qe. The charge difference is small enough to have a negligible effect on the trajectories, but large enough to be detected by the filter. Figure 1 shows the settings used to create the trajectory plot of Fig. 2. Note that the 10.0 keV electrons are shown in blue and the 12.0 keV electrons in red.

Trajectory plot with color-coding by injection energy

Figure 2. Trajectory plot with color-coding by injection energy.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.

Multiple particle species in Trak

Trak can handle multiple charged-particle species in a single run. The program can even determine self-consistent space-charge-limited emission for different particle types from multiple emission surfaces. When responding to a user question about run setups, I realized that there was no example to illustrate the capability in the Trak library. This example describes the example BIPOLAR that I added to the code distribution.

To demonstrate the validity and accuracy of a program, the best example is one with a known answer. In this case, I modeled bipolar flow, counter-streaming electrons and protons in a planar gap. Section 6.4 of my book Charged Particle Beams (http://www.fieldp.com/cpb.html) reviews the theory. It states that 1) the ratio of current densities is jp/je = √(me/mp) and 2) the current density of each species is 1.86 times the standard Child-law prediction because of mutual space-charge cancellation.

I set up a planar solution (variations in x-y, infinite length in z) with a gap d = 1.0 cm from x = 0 cm to x = 1.0 cm. The applied voltage was V0 = 10.0 kV. The solution volume extended from y = 0.0 cm to y = 1.0 cm with Neumann boundary conditions at the top and bottom. Regions 2 and 3 represented the physical cathode and anode and Regions 4 and 5 were emission line regions on their surfaces. There was one tricky point in creating a 1D solution — Trak does not accept emission points on the solution-volume boundaries. As a workaround, I added thin elements on the top and bottom and defined the emission surfaces so they extended only to the node just inside the boundary.

The other issue was to use the long form of the Dt command for efficient tracking of particles with a large difference in mass:

Dt  DtRef  Mass

The quantity DtRef is the time step for a particle with mass 1 AMU (the proton). The velocity of a proton at 10 keV is 1.38E6 m/s. The time to cross an element of width 0.0001 m was 7.24E-11 s. I used a value DtRef = 5.0E-11 s..

The Trak calculation gave a total current (per meter in z) of 433.2 A. Using the formulas in Charged Particle Beams, the Child limit for bare electron and proton beams is je = 2.33 A/cm2 and   jp = 0.0544 A/cm2. The expected total current per meter is therefore

1.86*(2.33+0.0544)*100.0 = 443.5 A/cm2,

within 2.2% of the code result. Using filtering by x position in GenDist to separate electrons and protons, I found the following values for the proton and electron current per meter: Ip = 10.119 A and Ie = 430.34 A. The ratio of currents was Ip/Ie = 0.0235, close to the theoretical value of 0.0233.

This zip file contains the input files for this example: bipolar.zip.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.