# Modeling ferromagnetic materials in a primitive python FEA

I am working on a primitive 2D DIY magnetic "FEA" solver for a class project. It's really just a big NumPy array where each element in the array can be one of the following:

• Air
• Current carrying wire with current into or out of the screen (Biot-Savart applied to every other element in the array)
• A permanent magnet (current sheets into and out of the screen on the sides of the magnetic boundary parallel to the magnetization vector, so again Biot Savart applied to every other element in the array but this time in sheets of current)

These are working well ath the moment and I can create nice vector plots showing $$\vec{B}$$ across the entire array.

Now I am trying to add functionality for ferromagnetic materials. My approach for this has been to run a loop that goes like this:

1. Ignore ferromagnetic regions and just calculate $$\vec{B}$$ everywhere in the array. Store this as B_first
2. Now in ferromagnetic regions, calculate the induced magnetization vector $$\vec{M}$$ equal to $$\vec{B}$$ * susceptibility / $$\mu_{\text{material}}$$
3. Calculate the $$\vec{B}$$ field created by this newly magnetized material (as though it were a permanent magnet) and add it to B_first
4. Create a vector plot showing this new field
5. Return to step 2, but now use the updated field to calculate a new induced magnetization vector...

From some talks with more senior folks in our group I was under the impression that this would converge and I would get a nice series of vector plots showing everything converging. As a test case I am trying to visualize a magnet with $$\vec{M}$$ parallel to an adjacent block of steel getting its flux "absorbed" by the steel yoke.

At the moment my solution is not converging and the field inside the steel is growing pretty linearly with successive iterations, and eventually becoming many orders of magnitude larger that the original permanent magnet.

Is anyone familiar with this sort of simulation and potentially aware of what could be going wrong and interested in taking a look? If so I can send code/plots and would be very keen to discuss more.

For ferromagnetic materials such as metallic iron, $$\vec{B}=\mu \vec{H}$$ should hold. From the relation $$\vec{B}=\mu \vec{H}$$, $$\vec{B}$$ and $$\vec{H}$$ are parallel inside the ferromagnet (remember that $$\vec{B}$$ and $$\vec{H}$$ are approximately anti-parallel inside magnet). This seems a somewhat delicate fact, and in my experience it is not easy to achieve this state of affairs within numerical calculations of integration type ( but probably, there may be people who have done it.).