Running a 3D electrostatic solution: mesh generation and solution

This article is part of the continuing online course, Electric and magnetic field calculations with finite-element methods. This article continues an example 3D electrostatic solution. The previous article introduced the concept of mesh Parts and Regions. The emphasis was on input from 3D CAD programs via STL files. Here, we’ll consider how MetaMesh uses the STL data and processes extra information to create the result of Fig. 1. We could go through all the interactive operations in Geometer to define the parts, but there’s a quicker way to understand the setup. The end result Geometer activity is a succinct MetaMesh input script that includes all the information. We’ll study the script directly. A following article shows how to use Geometer to create a solution from scratch.

Conformal mesh created by MetaMesh

Figure 1. Conformal mesh created by MetaMesh, outline of the VACUUM region.

Run AMaze and launch MetaMesh. We’ll assume that the Data folder is still pointing to the working directory of the previous article. Click File/Load MIN file and choose sheetgun.min. Then click File/Edit MIN file to view the contents of the script. The first part of the Global section has detailed specifications for element sizes over multiple zones along the axes. The feature was necessary for accuracy in subsequent beam emission calculations – the details are not of immediate concern. Instead, we’ll concentrate on the Part definitions shown in Fig. 2 [1]

MetaMesh input script sheetgun.min

Figure 2. MetaMesh input script sheetgun.min.

The Parts named Anode, Focus and Cathode (marked in green) use the STL models that we previously discussed. The Type command of each part section reads the associated STL file. The Fab (fabrication) command defines fitting controls — default values are usually sufficient. Although the intent is build the mesh around the assembly positions of the STL parts, we do want to reverse the emission direction to +z. This is accomplished with the commands

Rotate 180.0 0.0 0.0 XYZ

Each command rotates the coordinates of all facet nodes 180° about the x-axis relative to the assembly origin. The commands

Surface Region VACUUM

instruct the code to adjust all element facets of the part adjacent to Vacuum elements to conform closely to the defined surface of the STL model.

The assembly is located inside a cylindrical vacuum chamber of radius 1.23″ that extends from z = 0.0″ to the boundary of the solution volume at z = -1.50″. For high speed, HiPhi uses a structured conformal mesh, which means that the elements are arranged like a big apartment building. The term conformal means that apartments are not necessarily right-angle boxes [2]. The solution volume of a structured mesh is always a box. To create the vacuum chamber, we fill the entire volume with elements assigned to the fixed-potential region Ground (part VChamber) and then carve out a cylindrical volume of elements associated with the Vacuum region. The part ExitPort is a slot to accommodate the beam exit aperture. The final part (Support) is a cylindrical pipe attached to the cathode. Note that the part belongs to the same Region as the cathode and the focus electrode.  The physical property of the Region in the HiPhi solution is is a fixed potential (-74.0 kV).

Exit the text editor and click the Process mesh command. MetaMesh reports the progress of the mesh creation. When the program is finished, right-click to exit the report and use the File/Save mesh command to create the file sheetgun.mdf (mesh definition file).

HiPhi dialog

Figure 3. HiPhi dialog to set program parameters and the physical properties of regions.

We can now proceed to the HiPhi solution. Run the program, click Setup and choose sheetgun.mdf to open the dialog of Fig. 3. Fill in the values as shown. Click OK and save the results as sheetgun.hin. The resulting file has the contents:

Mesh = sheetgun
DUnit =   3.9370E+01
ResTarget =   5.0000E-08
MaxCycle =    4000
Parallel = 4
* Region 1: GROUND
Potential(1) =   0.0000E+00
* Region 2: VACUUM
Epsi(2) =   1.0000E+00
* Region 3: CATHODE
Potential(3) =  -7.4000E+04

The clear format of the file makes it easy to identify the functions of the commands. I have the Pro version of HiPhi and a core i7 computer, so I added the command

Parallel = 4

with a text editor to speed things up.

Click the Run button and choose the input file to generate a solution recorded as the binary file sheetgun.hou. We can get some useful information by checking the text listing files and sheetgun.hls created by MetaMesh and HiPhi. For example, I found that the mesh had 1,295,952 active [3] elements and took 19 seconds to process on my computer. The 4000 matrix-inversion cycles in HiPhi took 363 seconds to reach a relative residual of 8.6E-07 [4].

To complete the exercise, run PhiView and load sheetgun.hou. You can experiment with the options for 2D and 3D plots. Figure 4 shows an example. Here, I wanted to get a quick overview of levels of |E| on the focus-electrode surface. The red areas correspond to |E| = 16.7 MV/m (167 kV/cm).

3D plot of the assembly

Figure 4. 3D plot of the assembly showing levels of |E| on the surface of the focus electrode.


[1] The attractive text display is from the ConText editor using syntax highlighting definitions that we supply with the software.

[2] The general term for a six-sided solid where facets may not be at right angles is a hexahedron.

[3] The term active means that the elements are not part of a fixed-potential region.

[4] The relative residual, a measure of the accuracy of the iterative solution, should be much smaller than unity.

[5] Contact us :

[6] Field Precision home page:


Running a 3D electrostatic solution: STL input

This article is part of the continuing online course, Electric and magnetic field calculations with finite-element methods. As in the introduction to the 3D software, we’ll run through a complete field solution for a sophisticated gun designed at the Lawrence Livermore National Laboratory to generate sheet electron beams . We’ll use prepared input files — later articles will help you prepare you own inputs.

To start, run the AMaze program launcher and set the Data folder command to point to a working directory. Copy the following files from the MetaMesh example directory to the working directory:


The text file sheetbeam.min is a set of instructions to MetaMesh giving the shapes and assembly method of the gun parts. The instructions were built in the interactive environment of Geometer. The text file sheetbeam.hin gives HiPhi the material properties and other parameters needed to generate a finite-element solution. The stereolithography files cathode.stl, focus.stl and output.stl were supplied by a laboratory engineer who exported them from a Solidworks model. They define shapes that are too complex to create from a summation of simple solid models (e.g., box, cylinder,…).

As we saw in a previous article, the procedure for 2D shapes is relatively simple. Objects are created from a boundary outline of line and arc vectors. The outlines may be derived from or exported to DXF files. In contrast, shapes encountered in the 3D world may be vastly more complex. A flexible approach is essential to make the task manageable. The following principles underlie the operation of Geometer and MetaMesh. Their implications will become clear as we progress through the example:

  • Objects (like electrodes and dielectrics) are constructed from one or more Parts.
  • Physical properties (like potential and relative dielectric constant) as assigned to Regions. Each Part belongs to a Region.
  • An object may be constructed from several Parts associated with the same Region. For example, a grounded electrode could be composed of Parts that belong to a Region that has a fixed potential 0.0 V.
  • Parts are processed in the order in which they appear in the MetaMesh input file. The currently-processed Part overwrites any shared volumes with previously-processed Parts.

Parts are defined at a standard position and orientation. They may be moved or rotated to a final configuration in the solution space.

It’s clear that the idea of a Part is central to the mesh building process. There are two types:

  • Parametric models for simple shapes like cones or spheres. These models are built into Geometer and MetaMesh and do not require external data.
  • Arbitrary shapes created with 3D CAD programs like SolidWorks and exported as files in the STL format. These files must be loaded into Geometer or MetaMesh.

A mesh may combine both types of parts.

Let’s get started with the calculation. To begin, it’s useful to understand the information in STL files. Run Geometer and click on the STL viewer command. The program becomes a full-featured STL viewer that can display the shapes defined by individual STL files and their relationships in space. Click on File/Add model and choose focus.stl. Click on View control/Orthogal/perspective to enter the 3D mode. This is an interactive environment based on OpenGL. Move the mouse cursor to the sides and left click to change the view. You can also move the cursor toward the center and left or right-click to zoom in or out. The view looks like the right-hand side of Fig. 1. The part is a dish-shaped focusing electrode with a hole for the cathode.

Geometer STL viewer

Figure 1. Geometer STL viewer, showing solid and wireframe views of focus.stl.

To see the inherent data of the STL file, click Plot control/Model display. Uncheck the Solid box to create a view like the left-hand side of Fig. 1. The view shows the set of contiguous triangular facets that define the surface of the part. The STL file is simply a list of facets. The facets shown are typical of those created by 3D CAD programs. Although all facets are theoretically correct, they may vary considerably in scale. To determine if an element is inside  an STL shape, MetaMesh must analyze all facets, an intensive activity where parallel processing is a big advantage.

Return to a solid view of the surface and load the other two parts to get a view like that of Fig. 2. Per my request, the mechanical engineer exported the parts from a SolidWorks assembly. In this case, the coordinates of facets in the STL file are absolute with respect to the assembly space rather than relative to the part. Therefore, the facets not only define the shape of parts, but also their positions and orientations relative to each other[1].

 Full set of STL shapes for the assembly

Figure 2. Full set of STL shapes for the assembly.

The parts as loaded define the full volumes of the cathode and focus electrode, but only half of the output transport tube. This is not a problem because by symmetry it is necessary only to do a calculation in the first quadrant of the x-y plane. MetaMesh automatically clips over-sized parts. There is one problem that you can see by shifting to the Orthogonal view and clicking View control/Y normal axis. The Solidworks model has the beam moving in the -z direction. To match transport calculations, we want the beam to move in +z following the standard convention. The modification will be made during MetaMesh processing.

Exit Geometer and run MetaMesh. Click File/Load MIN file and choose sheetgun.min. Click the Process mesh command and wait for the program to complete it’s analysis. It takes about 20 seconds to generate a mesh with 1.3 million elements on a multi-core computer. Right-click to close the messages. Choose Plot3D to observe the display of Fig. 3. The plot (which shows the surfaces of elements of the vacuum region) gives a good output of all the electrodes. The parts defined by STL files are present, but much more has happened. The mesh has been created in the first quadrant with fine resolution of the cathode surface, the hexahedron elements conform closely to the theoretical shapes, the beam direction points in +z, the assembly is inside a cylindrical vacuum chamber and there is a support rod for the cathode. All will be explained in the next article where we take a close look at the MetaMesh script.

 Completed mesh

Figure 3. Completed mesh, outline of the vacuum region.


[1] (Note that we could move the parts to different positions by adding shift operations in Geometer.

[2] Contact us :

[3] Field Precision home page:

Adding 3D software

With a good background in 2D finite-element solutions, it’s time to advance to the more challenging and sophisticated task of full 3D electrostatic and magnetostatic field calculations. In this article (the twelth in the online course Electric and magnetic field calculations with finite-element methods), we’ll make two assumptions:

  • You purchased or requested a trial of HiPhi, our unitized package for electrostatic and high-voltage devices.
  • You have already installed and activated the Basic version of EStat (2D electrostatics).

In this case, installation is quite simple. Activation of the computer will carry over from the existing programs.

If you make a purchase or request a trial, we will send an E mail message with a link to download HiPhi:

Package: HiPhi (Basic)
User: bin16
Password: HklW29567$

Download and run the installer. It adds the 3D programs to create the directory structure shown in Fig. 1. There are now two sub-directories in the main program directory: tricomp (the existing 2D software) and amaze (the new 3D programs). Note that the 2D and 3D programs are not directly coupled because they use completely different types of meshes (triangles in a plane versus hexahedrons in a volume). Note further that if you do not have an existing installation, it will be necessary to follow the process described in a previous article.

Directory setup with 2D and 3D programs

Figure 1. Directory setup with 2D and 3D programs.

Let’s review the contents of the amaze directory:

  • amaze.exe: A program launcher to run the HiPhi programs and to organize data files.
  • dielectric_constants.html and physcons.pdf: Useful reference material for electric field calculations.
  • geometer.exe: An interactive 3D environment based on OpenGL for building assemblies from simple shapes or exported parts from SolidWorks and other CAD programs.
  • hiphi.exe: The main FEM solution program.
  • hiphi.pdf: The HiPhi Instruction manual that can be called using tools in hiphi.exe and phiview.exe.
  • metamesh.exe and metamesh.pdf: Automatic conformal mesh generator, program and instruction manual.
  • notify.exe and notify.wav:  A beep program and waveform to support the Task feature of amaze.exe.
  • phiview.exe: A flexible graphical post-processor for analyzing HiPhi solutions.
  • phiview_conductive.cfg, phiview_dielectric.cfg and phiview_force.cfg  Standard PhiView configurations for different types of solutions and analyses.

Following article will clarify the functions of the components and their relationship to each other.

Amaze program launcher

Figure 2. Amaze program launcher.

Click on the desktop shortcut to Amaze created by the installer to run the program controller shown in Fig. 2. It contains a number of buttons to run the 3D programs and perform other tasks. The buttons are organized into the following groups:


  • Tasks: Advanced function for power users to create and to run batch files for automatic control of complex calculations.
  • Mesh generation: As we saw in the previous chapter on 2D calculations, defining the shapes of physical objects and dividing the solution space into small pieces is an essential part of the finite-element process.
  • Solution: The main FEM solution program. In Fig. 2, only HiPhi has been installed.
  • Analysis: Explore the solutions, creating plots and numerical results.
  • Utilities: Prepare inputs and display outputs for special solution programs. HiPhi does not use any of the utilities.
  • Control: Organize your data.
  • Tools: Call up text editors, terminal windows and other programs that you use in your work.

There is also an Activate button in case you haven’t previously installed 2D software. The next article describes a quick calculation that demonstrates how everything fits together.


[1] Contact us :

[2] Field Precision home page:

Magnetostatic solutions: permanent magnets

As the final topic in two-dimensional magnetostatics, we’ll consider solutions with permanent magnets. To start, it’s useful to review how permanent magnets work. Introductory electromagnetic texts often don’t treat the topic, and the explanations in many specialized references are overly complex.

The electrons of many atoms carry a circulating current, like a small current loop. Such an atom has a magnetic moment , a vector pointing from the center of the loop normal to the current. The vector lies in the the direction of the magnetic flux density B created by the loop. The distinguishing feature of ferromagnetic materials is that the atoms prefer to orient their currents in the same direction (a quantum mechanical effect). A region where all atoms are aligned is called a domain. Figure 1:top shows the summation of aligned atomic currents in a domain. Currents on the inside cancel each other, so the net result is a surface current in the direction normal to the magnetic moment.

Domains and atomic currents

Figure 1. Top: alignment of atomic currents in a domain. Bottom: alignment of domains in unmagnetized material.

In the natural state of a ferromagnetic material, the orientation of domains is randomized so that there is no macroscopic field outside the material (Fig. 1:bottom). External fields would require energy to generate. Let’s say that we supplied that energy by placing the material inside a strong magnet coil. In this case, the domains line up and the current of all domains sum up as in Fig. 1:top. The domain currents cancel inside, but there is a net surface current on the object. If we turn off the coil, the domains of a soft magnetic material return to a random distribution. On the other hand, suppose we could physically lock the domains in position before we turned off the coil. In such a hard material, the object retains its surface current and can generate external fields. The energy for the field was supplied by the magnetizing coil and it was bound in the material. Such a object is called a permanent magnet.

The distinguishing feature of modern permanent-magnet materials like neodymium-iron and samarium-cobalt is that domain locking is extremely strong. The domains remain lined up, independent of external processes. This property makes it simple to model the materials. Let’s review some definitions and facts. The direction of the magnetic moments of the aligned atoms (and domains) is called the magnetization direction. Suppose we have a permanent magnet that is long along the direction of magnetization and self-connected (e.g. a large torus). In this case, there is no external field and the flux density inside the material is generated entirely by the surface currents. This intrinsic flux density is called the remanence flux, Br. A typical value for neodymium iron is Br = 1.6 tesla.

We can express the surface current density in terms of the remanence flux. Take Br as a vector pointing along the direction of magnetization and let ns be a vector normal to the surface of the permanent magnet. The surface current density is given by

Js = (Br × ns)/μ0 (A/m).  [1]

We can use PerMag to confirm this physical interpretation. Figure 2:top shows a calculation in cylindrical geometry for an annular permanent magnet in space (i.e., no iron, coils or other permanent magnets). The remenance field is Br = 1.5 tesla and the direction of magnetization is along z. The plot shows lines of B. According to Eq. 1, we would get the same results by replacing the permanent magnet with thin surface current layers. By the properties of the cross product in Eq. 1, there should be no current on the ends and uniform current density Js = 1.193 Ma/m on the inner and outer radial surfaces. In the second model (Fig. 2:bottom), the permanent magnet is replaced by two thin solenoid coils (length = 0.08 m, thickness = 0.00125 m). The total current of the outer coil is 1.193E6*0.08 = 95470.0 A and the current on the inner coil is -95470.0, Figure 2:bottom shows that the calculated lines of B are indistinguishable from the permanent magnet calculation. For a quantitative check, we can compare scans of Bz along the axis. Figure 3 shows the results. The small difference is a result of the finite thickness of the layers.

Lines of magnetic flux density for an annular permanent magnet (top) versus equivalent surface current layers (bottom).

Figure 2. Lines of magnetic flux density for an annular permanent magnet (top) versus equivalent surface current layers (bottom).

Scans of By(x,0) for an annular permanent magnet versus equivalent surface current layers.

Figure 3. Scans of By(x,0) for an annular permanent magnet (line) versus equivalent surface current layers (red symbols).

If permanent magnets are that simple, where could confusion arise? Unfortunately, things are more challenging when we talk about older materials like Alnico. Such materials are characterized by a demagnetization curve. The curve is usually a plot of B inside the permanent magnet versus H, where H is the applied magnetic field that arises from coils and other permanent magnets. It is easier to see the meaning of the curve if we make the plot in terms of the applied magnetic flux density, B0 = μ0*H. Applied fields have no effect on a modern materials. Therefore, the total flux density inside the material is simply

B = Br – B0.  [2]

Figure 4a shows a plot. The value of applied magnetic field where B = 0.0 is called the coercive force Hc. For modern materials, the coercive force is

Hc = -Br/μ0.  [3]

Demagnetization curves for modern (a) and older (b) permanent magnet materials.

Figure 4. Demagnetization curves for modern (a) and older (b) permanent magnet materials.

It is clear that neither the demagnetization curve or the coercive force are particularly meaningful for materials like NdFeB. On the other hand, the concepts are useful for older materials. Here, the domains are not tightly locked. Putting such a magnet in a device with an air gap causes degradation of the alignment so that the total flux density falls below the ideal curve (Fig. 4b). PerMag can solve such problems, but it is important to clarify what such solutions mean. Suppose you bought an Alnico magnet that was magnetized at the factory and shipped with an iron flux return clamp (i.e., B = Br). If you could move it instantaneously to the application device, then the operating point calculated by PerMag would apply. On the other hand, if you removed the clamp during assembly so the permanent magnet was exposed to an air gap and possibly dropped it on the floor, then the material would have undergone an irreversible degradation and the fields generated would be lower than the numerical prediction. The only way to get the theoretical performance based on the demagnetization curve is to magnetize the material in situ. One big advantage of modern materials (beside their strength) is that they achieve the predicted performance, even if they are magnetized at the factory, shipped from China and moved to the assembly.

We’ll conclude with an application calculation, a bending magnet for an ion spectrometer. This is a long assembly, and we will do a preliminary 2D calculation in a cross section. Figure 5:top shows half the geometry (there is a symmetry plane with a Neumann boundary at y = 0.0 cm). The goal is to create a dipole field By in the air gap to bend ions moving in z in the x-direction. The arrow shows the magnetization direction of the permanent magnet (NdFeB with Br = 1.6 tesla). Two features ensure the maximum air-gap field for the given magnet surface currents:

The figure shows the calculated lines of B. Figure 6 shows a scan of By along x across the air gap. The field is strong but non-uniform, probably of little use in a spectrometer. Here, where we can take advantage of one of the properties of soft steel, field shaping. Consider the effect of adding a steel layer adjacent to the gap. Figure 5:bottom shows the change in lines of B. The steel shifts flux away from the center out to the edges of the gap. Figure 6 shows the effect on the field scan. With some sacrifice in the magnitude of the flux, the steel insert gives a working volume of approximately uniform field.

Application example, permanent-magnet ion spectrometer

Figure 5. Application example, permanent-magnet ion spectrometer. Top: bare magnet. Bottom: magnet with steel insert for field shaping.

Spectrometer application, scan of By(x,0)

Figure 6. Spectrometer application, scan of By(x,0) with and without the steel insert.


[1] Contact us :

[2] Field Precision home page:


Magnetostatic solutions: when steel gets complicated

Material properties are the reason why magnetostatic solutions are generally more involved than electrostatic solutions. Dielectrics can usually be characterized by a single value of relative dielectric constant εr up to the point where they break down. On the other hand, magnetic materials do interesting things even at normal values of flux density B:

  • Magnet steels may exhibit variations of μr, with a drop to μr ≈ 1 at high field levels (saturation).
  • In permanent magnets, the domains are locked in place, almost independent of the applied field.

We’ll talk about how to model permanent magnets in the next article. Here, we’ll concentrate of non-linear magnetic materials.

First, some definitions. The quantity B0 is the magnetic flux density created by coils or permanent magnets. We’ll call B0 the applied flux density. In the absence of steel, the total flux density is B = B0. In the PerMag program, the relative magnetic permeability is defined as μr = B/B0. Therefore, μr = 1.0 in air or vacuum. In steel, the alignment of domains creates a material flux density that adds to the applied flux density, so that μr is greater than 1.0. As we saw in the previous article, the utility of steel follows from that fact that μr » 1.

Magnetic materials are characterized by curves of the form B versus B0 (or B versus H, where H = B0/μ0). If the variation is a straight line, then μr has a constant value and we say that the material is linear. Magnetic materials may exhibit linear behavior at low fields, but they always become non-linear at high fields (a few tesla). Equivalently, we can characterize materials with a plot of B versus μr. Figure 1 shows such a plot for Steel 1020, a common material for magnet cores. At low B, the value of μr is much greater than unity and varies considerably (i.e., the B-B0 curve is not a straight line). Physically, the curve implies that it takes some pushing to start aligning domains but things get easier when they begin to come around.

Plot of relative magnetic permeability versus B for Steel 1020

Figure 1. Plot of relative magnetic permeability versus B for Steel 1020.

A notable feature of the curve is that μr approaches 1.0 when all the domains are aligned. In this case, the material contribution to B does not increase as B0 gets larger, so that the slope B/B0 approaches unity. The flux density at which all domains are aligned is called the saturation flux density. For Steel 1020, the value is Bs = 1.75 tesla. PerMag can perform self-consistent calculations for non-linear materials using B-μr data. The program simultaneously adjusts the value of μr in elements based on the present value of B as it performs the iterative matrix solution of the finite-element equations.

Let’s begin a calculation. The first issue is whether the variations of μr below saturation will make a significant difference in the results. In other words, do we need to worry about getting the B-μr curve exactly right at low field values? To illustrate, we’ll use the example of the H magnet illustrated in Fig. 2. It has a planar rather than cylindrical geometry. Here, there are variations in x-y but the steel core and coils extend an infinite distance in z. Such a calculation is often a good approximation to the central section of a long dipole magnet for bending particle beams. Both sets of coils produce a upward-directed flux that is conducted by the steel to the air gap, the working volume.

Geometry and calculated lines of B in the cross section of an H magnet

Figure 2. Geometry and calculated lines of B in the cross section of an H magnet.

To start, we’ll address the question: how does the value of μr in the core affect the field distribution? We’ll use the command-line parameter method discussed in a previous article, assigning fixed μr values over the range 5.0 to 10,000 and comparing values |B| in the air gap. Before beginning, its useful to recognize that the solution of Fig 2 involves more work than is necessary. The solutions in each of the four quadrants are mirror images. Alternatively, we could seek a solution only in the first quadrant by applying appropriate symmetry conditions along the boundaries x = 0.0 and y = 0.0. Lines of B are normal to the boundary at y = 0.0. Here, we could apply the Neumann condition (the natural condition for a finite-element solution). Lines of B are parallel to the other three boundaries, implemented by the Dirichlet condition Az = 0.0.

Figure 3 shows a plot of By in the magnet air gap at the midplane (0.0,0.0) and edge (2.0,0.0) of the magnet as a function of μr. The relative magnetic permeability has a strong influence at low values but little effect for large values. We can understand what’s happening by inspecting Fig. 4, a plot of lines of B. The top illustration shows a solution with high μr. In this case, the low-reluctance core carries most of the flux. It flows up between the coils and returns across the air gap. A small fraction takes a short-cut around the outside of the outer coil (leakage flux). The reluctance of the air gap cause the flux to spread out to increase the cross-section area (fringing flux). Note that the lines of B entering and exiting the core are almost normal to the surface. In contrast, the lower illustration shows the case with low μr. In this case, the core has a high reluctance, increasing the leakage flux. A reduced fraction of the flux is transported to the air gap, hence the reduced value of By(0,0). Finally, we’ll consider why there is little change in the solution when μr changes by a factor of 10 between 1000 and 100000. The important quantity in finite-element magnetic field calculations is not μr, but γ = 1.0/μr. A value γ = 1.0 represents air, and a value γ ≈ 0.0 is characteristic of unsaturated steel or iron. Both 1/1000 and 1/10000 are close to zero, so the change in μr makes little change in the solution.

Variation of By as a function of relative magnetic permeability

Figure 3. Variation of By(0,0) [blue] and By(2,0) [red} in the H magnet as a function of relative magnetic permeability in the steel core.

Lines of B at high and low values of relative magnetic permeability

Figure 4. Lines of B at high and low values of relative magnetic permeability.

We’ll conclude with a full non-linear calculation for the H magnet. This is the script HMagnet.PIN that controls the PerMag calculation:

Mesh = HMagnet
Geometry = Rect
DUnit = 1.0000E+02
ResTarget = 1.0000E-09
MaxCycle = 8000
Avg = 0.05
Update = 20 500
* Region 1: AIR
  Mu(1) =   1.0000E+00
* Region 2: COREUPRT
  Mu(2,Table) = steel1020_permag.dat
* Region 3: COILUPRTIN
  Mu(3) =   1.0000E+00
* Low current
* Current(3) =   1000.0
* High current
  Current(3) =   10000.0
  Mu(4) =   1.0000E+00
* Low current
* Current(4) =  -1000.0
* High current
  Current(4) =  -10000.0
* Region 5: BOUNDARY
  VecPot(5) =   0.0000E+00

Region 2 is the steel core. Instead of assigning a fixed value, the code reads the table of information plotted in Fig. 1. The coils are set up for two calculations at high and low current. The values of drive current were chosen to give solutions where the core steel is well below and well above saturation.

Because two iterative processes are adjusted simultaneously, program parameters must be adjusted to ensure solution convergence. Some experimentation may be necessary. Consider the following commands:

MaxCycle = 8000
Avg = 0.05
Update = 20 500

The quantity MaxCycle is the maximum number of iterations, set to a fairly high value. The quantity Avg controls averaging of μr values between cycles. The low value is used to avoid oscillations. Finally, the Update command states that μr should be recalculated every 20 iteration cycles and that the program should wait 500 cycles before applying corrections to ensure a good starting solution.

Figure 5 shows the variation of μr at low and high current. At low current (upper solution), the relative permeability exceeds 1000 over most of the core volume. The low values at the top reflect the fact that |B| is small and the elements are on the left-hand side of the curve of Fig. 1. An inspection of lines of B in Fig. 6(top) shows that they are contained mainly in the core. Lines outside the core are normal to the surface. The midplane field is By(0,0) = 0.2459 tesla. At high current (Fig. 6 bottom), many regions of the core are driven into saturation, particularly the vertical piece adjacent to the air gap. The figure shows enhanced fringing flux near the air gap. The central field value is By(0.0) = 1.484, only 6.03 times the field at one tenth the drive current. Note that the field lines near the air gap are not normal to the core surfaces, affecting the profile of B across the gap.

Figure Spatial variation of relative magnetic permeability in the H magnet solution

Figure 5. Spatial variation of relative magnetic permeability in the H magnet solution for low (top) and high (bottom) drive currents.

Lines of B in the H magnet solution

Figure 6. Lines of B in the H magnet solution for low (top) and high (bottom) drive currents.


In summary, this article addressed the following topics:

  • Typical variations of μr for soft magnetic materials.
  • Physical implications of the magnitude of μr, particularly as affects the reluctance of steel structures in magnets.
  • Reducing the run time of a solution using symmetry boundaries.
  • How to set up and to interpret non-linear solutions in PerMag.

The next article will discuss how permanent magnets work and how to build solutions for 2D permanent magnet devices.


[1] Contact us :

[2] Field Precision home page: