Analysis of an Eddy Current Brake with FEMM

David Meeker
dmeeker@ieee.org

<<c:\\femm42\\mathfemm\\mathfemm.m

Preliminary Definitions

mm=0.001;
Inch=0.0254;
MS=10^6;
MA=10^6;
m=1;
(Subscript[μ, o])=4. Pi*10^(-7);
Hz=2*Pi;

Problem Parameters

A conductive plate moves within a double-sided array of steel-backed permanent magnets. This notebook computes the force of the brake vs. plate velocity.
One wavelength of the geometry, neglecting the end effects, looks like so:

BHmax = 40;   (* Energy product of the magnet  in MGOe *)
pp = 25*mm; (* Pole pitch of the magnets *)
wm = 20*mm; (* Magnet width *)
tm = 5*mm; (* Magnet thickness *)
tb  = 5*mm; (* Back iron thickness *)
tp = 2.5*mm; (* Conductive plate half-thickness *)
g = 1*mm; (* Air gap length between magnet and plate *)
σ = 30 ; (* Plate conductivity in MS/m *)
h = 100*mm; (* Depth of magnet array *)
ϵ = 20*mm; (* Plate overhang on either side of the magnets *)

Derived Quantities

Mechanical frequency of the magnet array

β=π/pp;

A big problem with using a 2D solver to analyze this sort of machine is that 3D transverse edge effects can significantly influence the performance of the brake.  Fringing near the transverse edge means that some of the PM flux doesn't link the plate, lessening the maximum force produced by the brake.  The resistance of the edge of the plate, where the currents "turn around", makes the plate look more resistive and pushes out the speed at which the peak force occurs.  Here, we will put in some kludgy corrections to account for these effects.

The total stack-up of unit permeability materials between plate centerline and iron is gtot:

gtot=(tp+g+tm);

Use this distance to get an appoximate derating to account for fringing at the transverse edges. This is an approximate kludge that is probably no good for an array that isn't wide compared to gtot. Say that the width of a fringing region is hf:

hf=gtot/2;

Shorten the active region by a bit, snipping off the fringing region of width hf on each side of the stator:

heff=h-2*hf;

Add the width of the fringing region to the width of the "end turn" region of the plate

ϵeff=ϵ+hf;

Correct the conductivity to account for the transverse edges of the plate

σeff = σ/(1+2/((β^2)*heff*ϵeff))

Remanence of the PMs in Tesla (see http : // www.femm.info/Archives/misc/BarMagnet.pdf)

Br=Sqrt[BHmax]/5.

Coercivity of the PMs in A/m

Hc=Br/(Subscript[μ, o])

You can replace a permanent magnet by an equivalent current sheet that produces the same external fields. For this particular case, the equivalent current sheet representation looks like:

We want to replace the current sheet by a sinusoidally distributed band of current density.  Most of the force in an eddy current brake is due to the fundamental, and if we just consider the fundamental, we can make a simulation in FEMM that includes "motion".

Find amplitude of the fundamental in a Fourier Series (http://en.wikipedia.org/wiki/Fourier_series) representation of the equivalent current density of the PMs in units of A/m^2.
There is a + current sheet of strenght Hc located at (pp-wm)/2 and a - current sheet located at  (pp+wm)/2.  The result is in units of A/m^2

Jmag=(Hc*Cos[β*(pp-wm)/2]- Hc*Cos[β*(pp+wm)/2])/Integrate[Cos[β x]^2,{x,0,pp}]

The analysis approach is to analyze the brake from the point of view of an observer attached to the plate.  From the point of view of this observer, the magnet's sinusoidally distributed current bands look like simple AC currents.  We can analyze the problem in FEMM by defining a distribution of currents in the region where the permanent magnets live that looks like Jmag*exp(I*β*x) where x indexes the position down the array.  The apparent frequency is ω=β*v, where v is the plate's velocity.

Draw Domain

This section uses the MathFEMM interface betwen Mathematica and FEMM to programmatically draw the brake geometry.

OpenFEMM[]

NewDocument[0]

Define Materials

MIGetMaterial["Air"];

MIAddMaterial["Steel",5000,5000];

MIAddMaterial["Plate",1,1,0,0,σeff ,0,0,1,0,0,0,1,1];

Draw Plate

MIDrawRectangle[{0,0},{2 pp,tp}];
MIAddBlockLabel[{pp,tp/2}];
MISelectLabel[{pp,tp/2}];
MISetBlockProp["Plate",1,0,"<None>",0,1,1];
MIClearSelected[];

Draw Air Gap

MIDrawRectangle[{0,tp},{2*pp,tp+g}];
MIAddBlockLabel[{pp,tp+g/2}];
MISelectLabel[{pp,tp+g/2}];
MISetBlockProp["Air",1,0,"<None>",0,1,1];
MIClearSelected[];

Draw Magnets

n=12;
For[k=1,k≤n,k++,
xl=2*(k-1)*pp/n;
xr=2*k*pp/n;
xm=(xl+xr)/2;
MIAddMaterial["J"<>ToString[k],1,1,0,(Jmag/MA)*Exp[I*β*xm],0 ,0,0,1,0,0,0,1,1];
MIDrawRectangle[{xl,tp+g+tm},{xr,tp+g}];
MIAddBlockLabel[{xm,tp+g+tm/2}];
MISelectLabel[{xm,tp+g+tm/2}];
MISetBlockProp["J"<>ToString[k],1,0,"<None>",0,0,1];
MIClearSelected[];
];

Draw Backiron

MIDrawRectangle[{0,tp+g+tm},{2*pp,tp+g+tm+tb}];
MIAddBlockLabel[{pp,tp+g+tm+tb/2}];
MISelectLabel[{pp,tp+g+tm+tb/2}];
MISetBlockProp["Steel",1,0,"<None>",0,0,1];
MIClearSelected[];

Add Boundary Conditions

MIAddBoundProp["A=0",0,0,0,0,0,0,0,0,0];
MIAddBoundProp["pbc1",0,0,0,0,0,0,0,0,4];
MIAddBoundProp["pbc2",0,0,0,0,0,0,0,0,4];
MIAddBoundProp["pbc3",0,0,0,0,0,0,0,0,4];
MIAddBoundProp["pbc4",0,0,0,0,0,0,0,0,4];

MISelectSegment[{pp,tp+g+tm+tb}];
MISetSegmentProp["A=0",0,1,0,0];
MIClearSelected[];

MISelectSegment[{0,tp+g+tm+tb/2}];
MISelectSegment[{2*pp,tp+g+tm+tb/2}];
MISetSegmentProp["pbc1",0,1,0,0];
MIClearSelected[];

MISelectSegment[{0,tp+g+tm/2}];
MISelectSegment[{2*pp,tp+g+tm/2}];
MISetSegmentProp["pbc2",0,1,0,0];
MIClearSelected[];

MISelectSegment[{0,tp+g/2}];
MISelectSegment[{2*pp,tp+g/2}];
MISetSegmentProp["pbc3",0,1,0,0];
MIClearSelected[];

MISelectSegment[{0,tp/2}];
MISelectSegment[{2*pp,tp/2}];
MISetSegmentProp["pbc4",0,1,0,0];
MIClearSelected[];

MIZoomNatural[]

Calculate Force vs Speed

Expected speed at which maximum braking force occurs in m/s, based on 1D linear induction motor theory:

ωc=((β^2)*gtot)/((σeff*MA/m)*(Subscript[μ, o])*tp)

vopt1d=ωc/β

We can use the quick analytical estimate of vopt to select the range over which we evaluate the performance of the brake.

MISaveAs[NotebookDirectory[]<>"temp.fem"];

data={};
For[v=0,v<2*vopt1d,v++,
MIProbDef[β*v/Hz,"meters","planar",10^(-8),heff,30];
MIAnalyze[];
MILoadSolution[];
MOGroupSelectBlock[1];
f=MOBlockIntegral[11];
data=Append[data,{v,f}]
];
v=.

CloseFEMM[]

ListPlot[data,Joined->True,ImageSize->500,Frame->True,GridLines->{Table[k,{k,0,25,5}],Table[-k,{k,0,600,100}]},FrameLabel->{"Speed, m/s","Force, N"},PlotStyle->{Directive[Thick]},BaseStyle->{FontFamily->"ComicSansMS",FontSize->14}]

opt=FindMinimum[Interpolation[data][vopt],{vopt,5,15}]

Fmax=2 opt[[1]];
vopt=vopt /. opt[[2]];

Theoretically, the braking force out to have the analytical form :

Fbrake=Fmax*(v/vopt)/(1+(v/vopt)^2)

Show[ListPlot[data,PlotRange->All],Plot[Fbrake,{v,0,2 vopt1d},PlotRange->All],ImageSize->500,Frame->True,GridLines->{Table[k,{k,0,25,5}],Table[-k,{k,0,600,100}]},FrameLabel->{"Speed, m/s","Force, N"},BaseStyle->{FontFamily->"ComicSansMS",FontSize->14}]

Conclusions

A basic analysis of a "long stator" eddy current brake has been presented in this notebook.  End effects have been neglected.  Long stator end effects could be modeled by simply modeling the entire length of the plate in the present style of analysis, adding a few extra wavelengths of stator on either side of the plate.  End effects make a force ripple at a frequency of 2βv, but for a reasonably long brake, the amplitude of the ripple is relatively small and the change to the DC force is negligible. See, for example,
http://esrdc.mit.edu/library/ESRDC_library/Kirtley_The_End _Effect _in _Short _S _ 1.pdf
http://www.dtic.mil/dtic/tr/fulltext/u2/a447307.pdf
http://ieeexplore.ieee.org/iel5/27/5680667/05613198.pdf

However, for "short stator" machines (e.g. an array of magnets riding down an aluminum fin), the end effect are a bit more signficant. A solver that includes motion is really needed to analyze the end effects on a short stator machine.  See, for example:
http://www.jeet.or.kr/ltkpsweb/pub/pubfpfile.aspx?ppseq=31
http://onlinelibrary.wiley.com/doi/10.1002/eej.4390980407/abstract