I had an idea for a project to control the position a ferrous object in 3D space to give the effect that the object is levitating and floating around in patterns. To achieve this, I will need some method to apply a force on the object to control its position and a method to measure its position. The structure that came to mind to achieve this was a Helmholtz coil. To apply a force onto the object, the coils of a 3D Helmholtz coil can be energized. To measured the position of the object, the inductance of the coils could be measured.
A Helmholtz coil is a device that can be used to generate a near uniform magnetic field within two coil. This is achieved by spacing the coils apart at a distance equal to the radius of the coils. The inductance of the coils should be a function of the position of the ferrous object.
To begin this project, I wanted to model a Helmholtz coil parametrically to tune demission's and gain further insight into the coils characteristics. To accomplish this, I used FEMM which is a free open-source software that can simulate low-frequency electromagnetic and electrostatic problems. This software can be controlled through MATLAB thus allowing for parametric control of the model.
Below is a simulation of a Helmholtz coil of radius 81.5mm with 1A passing through the coils of 10 windings. This size was chosen only because I already made a model using Fusion360 and 3D printed the design. Note that the magnetic flux lines are near vertical when passing near the middle of the coils. It's also noteworthy that the magnitude of the magnetic flux density is also near uniform in this region as shown by the same color.
Below is the MATLAB code used to parametrically generate this model. The code could be reduced significantly but I wanted to put effort towards eventually making a MATLAB based material library for FEMM and using this code to parametrically control any parameter. To run the code, the mfiles in the FEMM install folder must be added to the MATLAB path.
clc
clear all
close all
addpath('D:\Programs\FEMM\femm42\mfiles')
%% General Options;
Options.UseAxis = [true false false];
%% FEMM Options
Options.Frequency = 0; % Frequency of the circuit
Options.Units = 'millimeters'; % Units of the model
Options.ProblemType = 'planar'; % 2-D Planar problem
Options.Precision = 1e-8; % Precision required for the solver
Options.Depth = 1; % Depth into the page
Options.MinAngle = 30; % Min angle sent to the mesh generator
Options.Circuit = 'icoil'; % Name of the circuit
Options.CircuitCurrent = 1; % Current in the circuit
Options.CircuitType = 1; % Series connected circuit
%% FEMM Materials, Air
FEMM.Material.Air.MaterialName = 'air';
FEMM.Material.Air.muX = 1; % Relative permeability in the x- or r-direction
FEMM.Material.Air.muY = 1; % Relative permeability in the y- or z-direction
FEMM.Material.Air.Hc = 0; % Permanet magnet coercivity in Amps/Meter
FEMM.Material.Air.J = 0; % Applied current density
FEMM.Material.Air.Cduct = 0; % Electrical conductivity of the material in MS/m
FEMM.Material.Air.Lamd = 0; % Lamination thickness in millimeters
FEMM.Material.Air.Phihmax = 0; % Hysteresis lag angle in degress, used for nonlinear BH curves
FEMM.Material.Air.lamfill = 1; % Fraction of the volume occupied per lamination that is actually filled with iron, default 1 (completely filled)
FEMM.Material.Air.LamType = 0; % 0 → Not laminated / laminated in plane
FEMM.Material.Air.Phihx = 0; % Hysteresis lag in the x-direction in degress
FEMM.Material.Air.Phihy = 0; % Hysteresis lag in the y-direction in degress
FEMM.Material.Air.nstr = 0; % Number of strands in the wire build. 1 for Magnet Wire
FEMM.Material.Air.dwire = 0; % Diameter of each of the wires constituent strands in millimeters
FEMM.Material.Air.Data = num2cell([FEMM.Material.Air.muX,...
FEMM.Material.Air.muY,...
FEMM.Material.Air.Hc,...
FEMM.Material.Air.J,...
FEMM.Material.Air.Cduct,...
FEMM.Material.Air.Lamd,...
FEMM.Material.Air.Phihmax,...
FEMM.Material.Air.lamfill,...
FEMM.Material.Air.LamType,...
FEMM.Material.Air.Phihx,...
FEMM.Material.Air.Phihy,...
FEMM.Material.Air.nstr,...
FEMM.Material.Air.dwire]);
FEMM.Material.Air.AutoMesh = false;
FEMM.Material.Air.MeshSize = 5;
FEMM.Material.Air.MagDir = 0;
FEMM.Material.Air.Group = 0;
FEMM.Material.Air.Turns = 0;
%% FEMM Material Coil
FEMM.Material.Coil.MaterialName = 'coil';
FEMM.Material.Coil.muX = 1; % Relative permeability in the x- or r-direction
FEMM.Material.Coil.muY = 1; % Relative permeability in the y- or z-direction
FEMM.Material.Coil.Hc = 0; % Permanet magnet coercivity in Amps/Meter
FEMM.Material.Coil.J = 0; % Applied current density
FEMM.Material.Coil.Cduct = 58; % Electrical conductivity of the material in MS/m
FEMM.Material.Coil.Lamd = 0; % Lamination thickness in millimeters
FEMM.Material.Coil.Phihmax = 0; % Hysteresis lag angle in degress, used for nonlinear BH curves
FEMM.Material.Coil.lamfill = 1; % Fraction of the volume occupied per lamination that is actually filled with iron, default 1 (completely filled)
FEMM.Material.Coil.LamType = 3; % 3 → Magnet Wire
FEMM.Material.Coil.Phihx = 0; % Hysteresis lag in the x-direction in degress
FEMM.Material.Coil.Phihy = 0; % Hysteresis lag in the y-direction in degress
FEMM.Material.Coil.nstr = 1; % Number of strands in the wire build. 1 for Magnet Wire
FEMM.Material.Coil.dwire = 0.51054; % Diameter of each of the wires constituent strands in millimeters
FEMM.Material.Coil.Data = num2cell([FEMM.Material.Coil.muX,...
FEMM.Material.Coil.muY,...
FEMM.Material.Coil.Hc,...
FEMM.Material.Coil.J,...
FEMM.Material.Coil.Cduct,...
FEMM.Material.Coil.Lamd,...
FEMM.Material.Coil.Phihmax,...
FEMM.Material.Coil.lamfill,...
FEMM.Material.Coil.LamType,...
FEMM.Material.Coil.Phihx,...
FEMM.Material.Coil.Phihy,...
FEMM.Material.Coil.nstr,...
FEMM.Material.Coil.dwire]);
FEMM.Material.Coil.AutoMesh = false;
FEMM.Material.Coil.MeshSize = 5;
FEMM.Material.Coil.MagDir = 0;
FEMM.Material.Coil.Group = 0;
FEMM.Material.Coil.Turns = 10;
%% FEMM Boundrary
% for calling mi_makeABC()
% Note: calling mi_makeABC() without parameters is OK and will use defaults
FEMM.Boundrary.Shells = 7;
FEMM.Boundrary.Radius = 0;
FEMM.Boundrary.CenterX = 0;
FEMM.Boundrary.CenterY = 0;
FEMM.Boundrary.EdgeType = 0; % Dirichley out edge type
% for calling mi_addboundprop()
FEMM.Boundrary.Name = 'ABC';
FEMM.Boundrary.A0 = 0;
FEMM.Boundrary.A1 = 0;
FEMM.Boundrary.A2 = 0;
FEMM.Boundrary.Phi = 0;
FEMM.Boundrary.Mu = 0;
FEMM.Boundrary.Sig = 0;
FEMM.Boundrary.c0 = 0;
FEMM.Boundrary.c1 = 0;
FEMM.Boundrary.BdryFormat = 0;
%% Coil Specifications
Coil.Radius = [81.5 91.5 101.5];
Coil.WindingRadius = 1.5;
%% Open a new FEMM model and assign document properties
openfemm;
newdocument(0); % Create new magnetic problem
mi_probdef(Options.Frequency,Options.Units,Options.ProblemType,Options.Precision,Options.Depth,Options.MinAngle);
%% Create the circuits
mi_addcircprop(Options.Circuit,Options.CircuitCurrent,Options.CircuitType);
%% Add meterials
mi_addmaterial(FEMM.Material.Coil.MaterialName,FEMM.Material.Coil.Data{:}); % Add material to the material list
mi_addmaterial(FEMM.Material.Air.MaterialName,FEMM.Material.Air.Data{:}); % Add material to the material list
%% Draw the geometery of the problem
mi_drawrectangle(-Coil.Radius(1)-Coil.WindingRadius,Coil.Radius(1)/2-Coil.WindingRadius,-Coil.Radius(1)+Coil.WindingRadius,Coil.Radius(1)/2+Coil.WindingRadius); % Draw the top left coil
mi_drawrectangle(-Coil.Radius(1)-Coil.WindingRadius,-Coil.Radius(1)/2-Coil.WindingRadius,-Coil.Radius(1)+Coil.WindingRadius,-Coil.Radius(1)/2+Coil.WindingRadius); % Draw the bottom left coil
mi_drawrectangle(Coil.Radius(1)-Coil.WindingRadius,Coil.Radius(1)/2-Coil.WindingRadius,Coil.Radius(1)+Coil.WindingRadius,Coil.Radius(1)/2+Coil.WindingRadius); % Draw the top right coil
mi_drawrectangle(Coil.Radius(1)-Coil.WindingRadius,-Coil.Radius(1)/2-Coil.WindingRadius,Coil.Radius(1)+Coil.WindingRadius,-Coil.Radius(1)/2+Coil.WindingRadius); % Draw the bottom right coil
mi_zoomnatural();
%% Create block label for each element
mi_addblocklabel(0,0); % Block bable for air
mi_addblocklabel(-Coil.Radius(1),Coil.Radius(1)/2); % Block lable for top left coil
mi_addblocklabel(-Coil.Radius(1),-Coil.Radius(1)/2); % Block lable for bottom left coil
mi_addblocklabel(Coil.Radius(1),Coil.Radius(1)/2); % Block lable for top right coil
mi_addblocklabel(Coil.Radius(1),-Coil.Radius(1)/2); % Block lable for bottom right coil
%% Assign block label properties for each element
mi_selectlabel(0,0); % Block bable for air
mi_setblockprop(FEMM.Material.Air.MaterialName,FEMM.Material.Air.AutoMesh,FEMM.Material.Air.MeshSize,'<None>',FEMM.Material.Air.MagDir,FEMM.Material.Air.Group,FEMM.Material.Air.Turns);
mi_clearselected;
mi_selectlabel(-Coil.Radius(1),Coil.Radius(1)/2); % Block lable for top left coil
mi_setblockprop(FEMM.Material.Coil.MaterialName,FEMM.Material.Coil.AutoMesh,FEMM.Material.Coil.MeshSize,Options.Circuit,FEMM.Material.Coil.MagDir,FEMM.Material.Coil.Group,FEMM.Material.Coil.Turns);
mi_clearselected;
mi_selectlabel(-Coil.Radius(1),-Coil.Radius(1)/2); % Block lable for bottom left coil
mi_setblockprop(FEMM.Material.Coil.MaterialName,FEMM.Material.Coil.AutoMesh,FEMM.Material.Coil.MeshSize,Options.Circuit,FEMM.Material.Coil.MagDir,FEMM.Material.Coil.Group,FEMM.Material.Coil.Turns);
mi_clearselected;
mi_selectlabel(Coil.Radius(1),Coil.Radius(1)/2); % Block lable for top right coil
mi_setblockprop(FEMM.Material.Coil.MaterialName,FEMM.Material.Coil.AutoMesh,FEMM.Material.Coil.MeshSize,Options.Circuit,FEMM.Material.Coil.MagDir,FEMM.Material.Coil.Group,-FEMM.Material.Coil.Turns);
mi_clearselected;
mi_selectlabel(Coil.Radius(1),-Coil.Radius(1)/2); % Block lable for bottom right coil
mi_setblockprop(FEMM.Material.Coil.MaterialName,FEMM.Material.Coil.AutoMesh,FEMM.Material.Coil.MeshSize,Options.Circuit,FEMM.Material.Coil.MagDir,FEMM.Material.Coil.Group,-FEMM.Material.Coil.Turns);
mi_clearselected;
%% Create open boundrary layers and show the mesh
mi_makeABC();
%% Save the ducement
mi_saveas('D:\Projects\3D Prints\Helmholtz Coil\temp.fem');
%% Analyze the soluation and load the results
mi_createmesh();
mi_zoomnatural()
mi_analyze;
mi_loadsolution;
mo_showdensityplot(1,0,1e-3,1e-7,'mag')
mo_showcontourplot(50,0,3.5e-5,'real')
Measured.Raw = mo_getcircuitproperties(Options.Circuit);
Measured.Current = Measured.Raw(1);
Measured.VoltageDrop = Measured.Raw(2);
Measured.FluxLinkage = Measured.Raw(3);
Measured.FluxPerCurrent = Measured.FluxLinkage / Measured.Current; % This is the same as inductance when in the absense of a permnent magnet or other coils
%%
disp(['Measured Inductance: ', num2str(Measured.FluxPerCurrent), 'H'])
%% Close FEMM
pause()
closefemm;Below is a rendering of the 3D Helmholtz coil designed in Fusion 360. The design was done such that the entire assembly can be 3D printed. Each coil is assembled in 2 parts and the entire structure is press fit together.
The next steps of the project include
Simulating the change in inductance of the coils as a function of position
Design a circuit to energize the coils and measure the phase
UPDATE 2/27/2021
Continuing off the preceding work, the simulation was modified to add a small pure iron sphere in the model and the position of the sphere was swept in the x and y direction. At each position, the inductance of the coil was measured. In order to measure the inductance, the frequency of the circuits had to be non-zero. The plot below shows the measured inductance as a function of position. From this plot, it becomes apparent that the x-position of the sphere has a negligible impact on the measured inductance relative to the y-position. This means that to measure the position of the sphere in 3D, measurements on all three of the axes of the structure will be required.
The inductance as a function of position was measured at several different frequencies to determine if a relationship between change in inductance as a function of frequency is also present. The minimum, maximum, and delta of inductance was measured as follows.
10MHz: min = 4.3696e-07, max = 4.4144e-07, delta = 4.4832e-09 H
1MHz: min = 4.387e-07, max = 4.4319e-07, delta = 4.483e-09 H
100kHz: min = 4.4332e-07, max = 4.4780e-07, delta = 4.4769e-09 H
10kHz: min = 4.4501e-07, max = 4.4931e-07, delta = 4.2928e-09 H
1kHz: min = 4.4868e-07, max = 4.4949e-07, delta = 8.1876e-10 H
These results show that above 10kHz, the delta in inductance in the y-direction has a minimal increase. The delta is also only about 1% of the maximum over the ±30mm span tested. This small change will provide insight into to resolution required to measure the phase of the signal driving the coils.
Next, the frequency was set back to 0Hz and the stress tensor force of the sphere was measured at each position. The plot below shows a quiver plot of the measured forces in newtons. This shows that the force exerted on the sphere as a function of position has a non-linear relationship to position in a circular pattern. Note that the magnitude of the force acting on the sphere is strongest near the coils and weakest in the center of the structure.
Comments