-
David Preiss authoredDavid Preiss authored
Electromagnetic FEA Simulations
In python we can create rudimentary 2D and 3D FEA simulations of current carrying wires and permanent magnets in free space using only two equations:
For current carrying wires we use Biot-Savart:
Current carrying wires get boundary conditions, and in this case a single value representing the _Idl _term. It is assumed that all wires are pointing directly into or out of the page, so the cross product Idl x r_hat' term will always use an angle of 90 deg. Also note that we are calculating magnetic field strength H and not flux density B, despite the equation above, so no u_0 term.
For permanent magnets we can use the equations for point source magnetic dipoles:
Magnetic dipoles get boundary conditions and a magnetization vector m, which can be applied to a magnetic region at any angle and magnitude within the XY plane. The dirac delta term is applied only if calculating magnetic field strength inside a magnetic material.
The rest is surprisingly simple and can be implemented in about 100 lines of code in Python. We can make a numpy array that computes the effect of any dipoles or current carrying wires on every location of a uniformly spaced grid of "elements." Each element is assigned a field strength in X and Y, and a material property such as air, magnet, or wire. By iterating across every field-producing element's effect on every other element, we can create vector fields representing magnetic field strength H like the one below, which shows a permanent magnet (red bounding box) with m pointing in -X, and a current carrying wire (green bounding box) with current flowing out of the page.
Examining the results above, there are some inaccuracies and limitations worth noting for improvement:
-
First let's note that there appears to be divergence on the left and right bounds of the permanent magnet where the field changes direction. I initially assumed that this was incorrect, but it turns out H fields are indeed divergent at the boundary of a magnet, and Gauss's law for magnetism applies to magnetic flux density B, but not field strength H. This should really be apparent from the fact that H = B/u_0 - M, where the magnetization vector M is just a uniform vector field that sources from one side of the magnet and sinks into the other, as shown in the image linked to above.
-
The field strength inside the permanent magnet appears to be several orders of magnitude stronger than outside the magnet, and the field strength vectors are all extremely parallel (this is the result of that dirac delta term noted earlier). This is actually incorrect as noted by Erik, and the reason is that the magnetic dipole equation is not valid for calculating the field strength contribution of an element to itself (r must be significantly greater than the radius of the Bohr magneton or else the dipole equation suggests a field strength approaching infinity as we divide by r=0). With that said, the above simulation should be valid for field strength outside of permanent magnets
-
Solving for systems with ferritic or "magnetizable" materials requires iteration and annealing, and therefore this simulation tool is limited to materials with permeability of free space. A magnetized ferromagnetic material proceeds to produce an H field of its own, which in turn changes the H field generated on the first solve from any wires or permanent magnets. This process must then be iterated through until the solution converges to some acceptable degree.
-
Finally one additional thing worth addressing is that we have not yet seen a reason for commercial EM FEA modeling systems to require an infinite "air domain bound" to guarantee accuracy. This suggests that their underlying physics engine is doing something completely different, and it would be interesting to try and address that moving forward.
Attempt #2 - Magnets as Current Sheets
Lots of changes to make after several very productive white-boarding sessions with Erik, who came up with the solution to the dipole problem outlined below. First is that in order to address problem #2 above, and to properly calculate the field within a magnetized material, we need to switch from the dipole equation to something that will properly calculate a magnetic material's self-contribution. To do this, we can think of making the switch to effectively modeling a magnet with Biot-Savart, and rather than using the "free current" associated with charge flowing through a wire, we will use "bound current" associated with electron orbits. We can relate this bound current to the magnetization vector we are familiar with by taking its curl.
Assuming that the magnetization vector is uniform in our square cross-section of material, we know that the curl will be zero everywhere except along the perimeters (where M changes to 0 outside the boundary), and therefore our bound current can be modeled as sheets of currents flowing along the perimeter of the element, into and out of the screen. We can then plug this bound current into the equation for the magnetic vector potential:
With the end goal of calculating a magnetic field, we can use the Helmholtz Decomposition, to write any vector field F as the sum of the gradient of its scalar potential, Phi (for which curl = 0), and the curl of its vector potential, A (for which divergence = 0).
Because Gauss' law states that B is free of divergence, we can represent the vector field B created by some magnetization vector M as the entirety of the B field - for this reason we want to go straight to B rather than H.
Finally if we want to solve for the B field from a current sheet, we can integrate over the length of the sheet to get the x and y components of the B fields at an arbitrary point [x,y], for a current sheet extending from y' = a to y' = b.
This integral is for a vertical current sheet located at x = 0, so if we want to move it around we can substitute in the difference for x, and to rotate it 90 degrees we would substitute the vertical x values for the horizontal y values and vise versa, and then do the same for the resulting B fields. If we want to solve for a magnet, we must place two of these sheets on either side of our magnetic material, and this will result in the signature toroidal dipole-field shape. We can then plug this equation into the simulation to solve for either entire regions of magnetic material at once (which vastly reduces the number of calculations required), as well as how I was previously solving calculating each element individually. The figure below shows the same permanent magnet modeled as a sum of unit elements and modeled as a single multi-element region, where the fields are identical for both.
During material assignment, permanent magnets are assigned a bulk area with identical M vectors, however ferromagnetic materials that subsequently become magnetized must be modeled on an element by element basis, due to the fact that their magnetization vectors are a function of their applied B field, and therefore non-uniform.
Returning to Earlier Simulation
Now if we model the same magnet and wire we tried with our previous dipole equation, we appear to get a much more reasonable result. Note that now we are modeling the B field as opposed to H so there is no divergence within the permanent magnet. Also note that aside from scaling, the results should be identical outside of the permanent magnet region
With all of that out of the way, if you now imagine the B field of a magnet, it becomes helpful to visualize current sheets rather than dipole moments. Another satisfying outcome of this result is that if you've ever used a tesla-meter and taken point source measurements of a cylindrical magnet, you will notice that the field is significantly stronger at the edges in comparison to the middle, which tracks well with these equations. See the plot below of a 1x25 permanent magnet, from my simulations, and the femm simulation of two thin disc magnets in green adjacent to it.