Newer
Older
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](https://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law):
![alt_text](images/BiotSavart.png "BiotSavart")
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](https://en.wikipedia.org/wiki/Magnetic_moment):
![alt_text](images/H_dipole.png "H_dipole")
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](https://en.wikipedia.org/wiki/Magnetic_moment#Internal_magnetic_field_of_a_dipole).
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:
1) 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](https://en.wikipedia.org/wiki/Magnetic_field#H-field_and_magnetic_materials) at the boundary of a magnet, and [Gauss's law for magnetism](https://en.wikipedia.org/wiki/Gauss%27s_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.
2) 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
3) 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.
4) 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"](https://www.comsol.com/blogs/course-modeling-electromagnetic-coils-in-comsol/) 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.
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.
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
![jBound](images/jBound.png "jBound")
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](https://en.wikipedia.org/wiki/Magnetic_vector_potential):
![Afield](images/Afield.png "Afield")
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).
![helmholtz](images/helmholtz.png "helmholtz")
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.
![AtoB](images/AtoB.png "AtoB")
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.
![BfieldFinal](images/BfieldFinal.png "BfieldFinal")
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.
![continuousDiscrete](images/continuousDiscrete.png "continuousDiscrete")
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
![DIYsim2](images/DIYsim2.png "DIYsim2")
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](https://www.kjmagnetics.com/blog.asp?p=magnets-are-weird) adjacent to it.
![longboi](images/longboi.png "longboi")
## Validation with Experimental Data
Because our model still has some limitations (we haven't accounted for three dimensionality / no real units applied / not sure what an accurate M would be etc.) it's really only valid to look at these simulation results from a relative perspective and not an absolute one. This fits well with another dataset that I already have of B_y field measurements at increasing distances from a planar coil at 1A and 2A. To make a roughly analogous model, I simulate 5 flat wires conducting current into and out of the page as shown in the images below.
![validationExperiment](images/validationExperiment.png "validationExperiment")
Below we plot B_y measured along the dotted white line in the top left image, and compare it to B_y simulated along the dotted black arrow in the top right image. We observe that they decrease similarly with increasing distance, and both roughly linearly proportional to current. These results are promising, for the finite method and Biot-Savart physics, but don't yet tell us anything about our PM and ferritic physics.
![validation](images/validation.png "validation")
## Modeling Ferritic Materials
For designing motors (assuming we stick with magnetostatics), the final piece remaining here is support for magnetizable/ferritic materials, so let's give that a shot now. If we first designate a region as being ferritic, each element in that region gets a new permeability which we will use to scale the magnetization vector induced by any applied B fields, such that M = B/u_material. We can then resolve for the global B field with this new induced magnet, and repeat this process until the solution converges. There's definitely something wrong with my code here however, as the field inside of the ferritic material alternates in vertical stripes after appearing to converge over 8 iterations.
![ferritic](images/ferritic.gif "ferritic")
So in the end more work to do here. We expect to see most of the permanent magnet's flux to recirculate through the steel, which would also reduce the flux outside of both materials (aside from the path into and out of the steel).
## More Things
### Running the Code
All of the code is stored in this repository. From the command line you can call pitFinalFEA.py, which plots the B vector field contributed to an elements by a magnetic material at some distance away.
The relevant physics functions are imported from pitFinalPhysics.py, where calc_magnet is used for both permanent magnets and magnetized ferritic materials.
```python
calc_biot_B(Idl, r) # returns Bx and By for a a current*length Idl at some distance and orientation r
calc_magnet_B(x_pos, y_pos, mx, my, hx, hy, ox, oy) # returns Bx and By for some magnetized material with magnetization vector [mx,my]
```
To create materials we import pitFinalMaterials.py where we can append lines of code with the origin, height, width, and magnetic properties of an element or region of elements.
```python
wireMaker(Idl,wire_w,wire_h,[x_center, y_center]) # Idl sign will dictate right hand rule
magnetMaker([mx,my], mag_w, mag_h, [mag_ox, mag_oy])
ferriteMaker(m,ferr_w1,ferr_h,[x_pos, y_pos]) # note that this will then create w*h 1x1 elements rather than monolithic block
```
In the validation section we talked about the limitations of looking at the simulation results in an absolute sense. If we wanted to, we should be able to extrapolate some symmetric 3D cases from these 2D results. Comsol has a very nice interface for this where you can actually revolve a 2D simulation to integrate a 3rd dimension without much computational overhead. I think
Where I think these simulations can get you into trouble is if you are modeling materials with different length scales. In that case the sum of contributions from each element dl all contribute to the field strength at point A, but with increasing effective radius, and decreasing angle theta, which results in the cross product dl x r_hat also decreasing in magnitude. From one of the psets I was under the impression that if we integrate the contribution of an infinite wire on a single point A, we scale Biot-Savart by a factor of 2r. From my understanding, if we assume all magnetic materials scale infinitely into the page, this would result in linear scaling and could still provide useful relative approximations.
In the case of a revolved coil for example, the center axis of the coil is at a fixed distance from the circumference of wire (assuming it isn't spiraling), and therefore we can scale Biot-Savart by 2*pi*r, but as we move towards the perimeter of the coil radially, the 3 dimensionality of the coil becomes more complex and nonlinear. I believe that unlike the infinite wire approximation, this would not result in linear scaling and perhaps would result in error unless accounted for?
We have a path to calculating forces via the [Lorentz Force Law](https://en.wikipedia.org/wiki/Lorentz_force) but I didn't get around to this.
* Software this project is trying to emulate:
* [FEMM - Finite Element Method Magnetics](https://www.femm.info/wiki/HomePage)
* [Comsol](https://www.comsol.com/support/learning-center/article/Modeling-Electromagnetic-Coils-8251)
* Erik pointed out that rather than go stright to B fields, most commercial solvers will likely calculate H and then... This explains why commercial solvers are so concerned with enforcing an infinite air boundary around their models...
* Real FEA packages do a lot of front-end work to mesh their models (triangular meshes rather than just enforcing a uniform grid). This allows for users to get higher resolution results in regions of interest, and avoid enforcing an element size dictated by the simulations smallest regions of interest. With that said, the meshing problem can be looked at as completely seperate from the physics, and I'd imagine at Comsol there is a meshing department and then various physics departments that don't interact too much with one another...
* Run on the GPU
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
# Motor Questions
The underlying motivation for this project was to answer some higher-order complexity questions pertaining to motor design. This section was written at the start of the project and uses mostly parametric sweeps from commercial FEA software.
![alt_text](images/FEAflux.png "image_tooltip")
## 1) What are the tradeoffs of using larger rotor magnets which provides more H, but increases the reluctance due to magnetic material's low permeability?
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
* Model an "infinitely long" 5mm diameter cylinder magnet and show diminishing returns for flux at the front of the magnet.
* Then model a coil with a yoke opposing it and increase the distance of the yoke to show diminishing returns there.
* Then with these two plots it should be possible to find an optimal distance, assuming nothing saturates. Some blogs have suggested that an infinitely long magnet is ideal. There must be flux contributed from both the rotor and stator to produce torque, but what is the right balance between these two? The plot below shows increasing magnet thickness (x-axis) and diminishing returns in shear force.
![alt_text](images/shear_v_magHeight.png "shearVmagHeight")
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
This plot jumps right to the conclusion and shows shear force with increasing magnet thickness along the x-axis. Starting at about 2 mm, we begin to see very diminishing returns in terms of shear force vs. increasing magnet thickness. To be determined is the tradeoff between the PM's contribution and the relative increase from
## 2) Is there a penalty to having a large insulator (which creates a high reluctance gap) between your motor steel and winding? Is it optimal to ensure that your winding is wrapped tightly against the inductor?
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
* H field produced is defined as the product of amps flowing through a winding and number of turns in that winding.
* Having a larger diameter coil increases the resistance of the winding, so perhaps the tradeoff is that we are creating a lot of flux that doesn't get used or gets canceled because there are more windings in the way.
* If you have a small steel core at the center of a large diameter of windings are you being punished?
## 3) What is the effect of filling an air core winding to the core with windings? Diminishing returns? Higher peak flux at center but lower exterior? Would an optimal coil have the smallest possible wall thickness?
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
* This would be a good one for FEA generally, but there's difficulty in that the ID increasing must effectively mean fewer coils, so it's necessary to compensate for both.
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
To start a simulation was run with constant amp-turns and increasing ID. Keep in mind that this was run with the optimal displacement as found in the sweep below for a 1mm ID at 2mm displacement. So it's possible that the magnet's location is convoluting this data a bit (more on that later). Interestingly there's a clear inflection point at 3.5mm corresponding to nearly double the shear force created in the .5mm case.
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
These results provide a strong basis for the possibility that tightly wound inner windings are contributing negatively, and so future sweeps here will need to be run manually to vary ID and have a % reduction in amp-turns corresponding with the coil area lost by increasing the coil's ID.
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
With that said, 3.5mm corresponds with the ID lying almost exactly on the center of the 3mm diameter of the magnet (see image below), so perhaps the takeaway is just that optimal shear will be produced at the point where the ID begins? To test this, lets re-run that offset code to make sure that the "Shear Displacement" parameter is still correctly set to 2mm with a now much larger ID.
![alt_text](images/annotated.png "image_tooltip")
![David Preiss David Preiss's avatar](/uploads/-/system/user/avatar/866/avatar.png?width=72)
David Preiss
committed
Again a sweep for optimal torque producing rotor displacement was run, analogous to load angle (max when 90 deg but unclear for this single pole linear displacement case). See CAD above for exact dimensions. Interestingly with this new ID, the optimal torque producing location has shifted up to 3.5mm. This suggests that iteration after tweaking certain parameters will be necessary to really converge on an optimal solution.
![alt_text](images/shear_v_displacement.png "shear_v_displacement")