%% Coreless wind turbine generator by Jose B. Almeida, February 2015 %% This OCTAVE script simulates an alternator with double magnet layers and %% coreless stator. %----------------------------------------------------------------------------------------------- %% modify these values to suit your needs units = 'millimeters'; % should be one of 'inches','millimeters', %% 'centimeters','mils','meters','micrometers' forceUnits = 'Newtons'; % should be one of 'lbf', 'Newtons', 'kgf' numberpp = 4; % number of pole pairs outerDia = 145.0; % external diameter of rotor magnetLength = 40.0; magnetThick = 10.0; magnetSep = 14.0; % separation between magnet layers magnetType= 'N52'; % has to be the EXACT name of a material in the %% materials library steelThick = 6; % thickness of magnetic circuit steel steelType = '1020 Steel'; % has to be the EXACT name of a material in the steelType = '1020 Steel'; % has to be the EXACT name of a material in the %% materials library bobbinReduce = 0.5; % Reduce the number of bobbins wireType = '1.79mm'; % has to be the EXACT name of a material in the %% materials library or add material wireDia = 1.79; % needs to be specified seperately %----------------------------------------------------------------------------------------------- addpath('C:/femm42/mfiles'); openfemm % opens an OctaveFEMM session %% We need to create a new Magnetostatics document to work on. %% Define the problem type. Magnetostatic; Units of mm; Axisymmetric; %% Precision of 10^(-8) for the linear solver; a placeholder of 0 for %% the depth dimension, and an angle constraint of 30 degrees newdocument(0); % the 0 specifies a magnetics problem mi_hidegrid(); %% define the problem type, equivalent to the top level menu Problem mi_probdef(0, units, 'planar', 1E-8, magnetLength); %% adds these materials from the Material Library to the project mi_getmaterial(magnetType); mi_getmaterial(steelType); % mi_getmaterial(wireType); mi_addmaterial(wireType,1,1,0,0,58,0,0,1,0,0,0,1,1.79); mi_getmaterial('Air'); %% add circuit properties for the 3 phases mi_addcircprop('A', 0, 1); mi_addcircprop('B', 0, 1); mi_addcircprop('C', 0, 1); %% this section is a bit of a pain. Although we have set the problem units above, %% that has no effect on the calculation of c0 for the %% Asymptotic Boundary Condition We therefore need to create a scaling factor %% for c0, dependant on the units c0_scale=1.0; if units=='micrometers' c0_scale= 1000000.0; elseif units=='millimeters' c0_scale= 1000.0; elseif units=='centimeters' c0_scale= 100.0; elseif units=='meters' c0_scale= 1.0; elseif units=='inches' c0_scale= 1.0/0.0254; elseif units=='mils' c0_scale= 1000.0/0.0254; end; forceScale = 1.0; if forceUnits=='Newtons' forceScale= 1.0; %% the natural units of the %% weighted stress tensor elseif forceUnits=='lbf' forceScale= 0.2248089; elseif forceUnits=='kgf' forceScale= 0.1019716; end; %% define the boundary uo = 4 * pi * 10-7; %% value of magnetic constant radius = outerDia * 0.7; %% set the window to a nice size for the problem mi_zoom(radius*1.2, 0.0, radius*1.5, radius*0.75); c0= c0_scale/(uo*radius); %% create the Asymptotic Boundary Condition for the problem mi_addboundprop('ABC',0,0,0,0,0,0,c0,0,2); mi_addnode(2*radius,-radius); mi_addnode(2*radius, radius); mi_addarc(2*radius,-radius,2*radius,radius,180,1); mi_selectarcsegment(2*radius,radius); %% make sure we set the Asymptotic boundary condition for the problem mi_setarcsegmentprop(1,'ABC',0,0); mi_addarc(2*radius,radius,2*radius,-radius,180,1); mi_selectarcsegment(-2*radius,radius); %% make sure we set the Asymptotic boundary condition for the problem mi_setarcsegmentprop(1,'ABC',0,0); %% draw outer and inner rotor shells mi_addnode(2*radius,-outerDia/2); mi_addnode(2*radius,outerDia/2); mi_addarc(2*radius,-outerDia/2,2*radius,outerDia/2,180,1); mi_addarc(2*radius,outerDia/2,2*radius,-outerDia/2,180,1); mi_addnode(2*radius,-(outerDia/2 - steelThick)); mi_addnode(2*radius,(outerDia/2 - steelThick)); mi_addarc(2*radius,-(outerDia/2 - steelThick),2*radius,(outerDia/2 - ... steelThick),180,1); mi_addarc(2*radius,(outerDia/2 - steelThick),2*radius,-(outerDia/2 - ... steelThick),180,1); innerDia = outerDia - 4*steelThick - 4*magnetThick - 2*magnetSep; mi_addnode(2*radius,-innerDia/2); mi_addnode(2*radius,innerDia/2); mi_addarc(2*radius,-innerDia/2,2*radius,innerDia/2,180,1); mi_addarc(2*radius,innerDia/2,2*radius,-innerDia/2,180,1); mi_addnode(2*radius,-(innerDia/2 + steelThick)); mi_addnode(2*radius,(innerDia/2 + steelThick)); mi_addarc(2*radius,-(innerDia/2 + steelThick),2*radius,(innerDia/2 + ... steelThick),180,1); mi_addarc(2*radius,(innerDia/2 + steelThick),2*radius,-(innerDia/2 + ... steelThick),180,1); mi_clearselected; mi_selectarcsegment(2*radius+outerDia/2,0); mi_selectarcsegment(2*radius-outerDia/2,0); mi_selectarcsegment(2*radius+innerDia/2,0); mi_selectarcsegment(2*radius-innerDia/2,0); mi_selectarcsegment(2*radius+(outerDia/2 - steelThick),0); mi_selectarcsegment(2*radius-(outerDia/2 - steelThick),0); mi_selectarcsegment(2*radius+(innerDia/2 + steelThick),0); mi_selectarcsegment(2*radius-(innerDia/2 + steelThick),0); mi_setarcsegmentprop(2.5, '', 0, 1); %% set outer cylinder properties mi_clearselected; %% Draw outer magnets magnetOuterArc = 180 / numberpp / 1.8; for i=0 : 2*numberpp - 1; mi_drawarc(2*radius - (outerDia/2 - steelThick - magnetThick) * ... cos(magnetOuterArc * pi/360 + i*pi/numberpp),(outerDia/2 - steelThick - ... magnetThick) * sin(magnetOuterArc * pi/360 + i*pi/numberpp),2*radius - ... (outerDia/2 - steelThick - magnetThick) * cos(-magnetOuterArc * pi/360 + ... i*pi/numberpp),(outerDia/2 - steelThick - magnetThick) * ... sin(-magnetOuterArc * pi/360 + i*pi/numberpp),magnetOuterArc,1); mi_drawline(2*radius - (outerDia/2 - steelThick - magnetThick) * ... cos(magnetOuterArc * pi/360 + i*pi/numberpp),(outerDia/2 - steelThick - ... magnetThick) * sin(magnetOuterArc * pi/360 + i*pi/numberpp),2*radius - ... (outerDia/2 - steelThick) * cos(magnetOuterArc * pi/360 + i*pi/numberpp)... ,(outerDia/2 - steelThick) * sin(magnetOuterArc * pi/360 + i*pi/numberpp)); mi_drawline(2*radius - (outerDia/2 - steelThick - magnetThick) * ... cos(-magnetOuterArc * pi/360 + i*pi/numberpp),(outerDia/2 - steelThick - ... magnetThick) * sin(-magnetOuterArc * pi/360 + i*pi/numberpp),2*radius - ... (outerDia/2 - steelThick) * cos(-magnetOuterArc * pi/360 + i*pi/numberpp)... ,(outerDia/2 - steelThick) * sin(-magnetOuterArc * pi/360 + i*pi/numberpp)); mi_selectarcsegment(2*radius - (outerDia/2 - steelThick - magnetThick) * ... cos(magnetOuterArc * pi/360 + i*pi/numberpp), (outerDia/2 - steelThick - ... magnetThick) * sin(magnetOuterArc * pi/360 + i*pi/numberpp)); mi_setarcsegmentprop(2.5, '', 0, 1); mi_selectsegment(2*radius - (outerDia/2 - steelThick - magnetThick) * ... cos(magnetOuterArc * pi/360 + i*pi/numberpp), (outerDia/2 - steelThick - ... magnetThick) * sin(magnetOuterArc * pi/360 + i*pi/numberpp)); mi_selectsegment(2*radius - (outerDia/2 - steelThick - magnetThick) * ... cos(-magnetOuterArc * pi/360 + i*pi/numberpp), (outerDia/2 - steelThick - ... magnetThick) * sin(-magnetOuterArc * pi/360 + i*pi/numberpp)); mi_setsegmentprop('', 0.5, 0, 0, 1); end %% Draw inner magnets magnetInnerArc = magnetOuterArc * (outerDia - 2*(steelThick+magnetThick))... /(innerDia + 2*(steelThick+magnetThick)); for i=0 : 2*numberpp - 1; mi_drawarc(2*radius - (innerDia/2 + steelThick + magnetThick) * ... cos(magnetInnerArc * pi/360 + i*pi/numberpp),(innerDia/2 + steelThick + ... magnetThick) * sin(magnetInnerArc * pi/360 + i*pi/numberpp),2*radius - ... (innerDia/2 + steelThick + magnetThick) * cos(-magnetInnerArc * pi/360 + ... i*pi/numberpp),(innerDia/2 + steelThick + magnetThick) * ... sin(-magnetInnerArc * pi/360 + i*pi/numberpp), magnetInnerArc,1); mi_drawline(2*radius - (innerDia/2 + steelThick + magnetThick) * ... cos(magnetInnerArc * pi/360 + i*pi/numberpp),(innerDia/2 + steelThick + ... magnetThick) * sin(magnetInnerArc * pi/360 + i*pi/numberpp),2*radius - ... (innerDia/2 + steelThick) * cos(magnetInnerArc * pi/360 + i*pi/numberpp)... ,(innerDia/2 + steelThick) * sin(magnetInnerArc * pi/360 + i*pi/numberpp)); mi_drawline(2*radius - (innerDia/2 + steelThick + magnetThick) * ... cos(-magnetInnerArc * pi/360 + i*pi/numberpp),(innerDia/2 + steelThick + ... magnetThick) * sin(-magnetInnerArc * pi/360 + i*pi/numberpp),2*radius - ... (innerDia/2 + steelThick) * cos(-magnetInnerArc * pi/360 + i*pi/numberpp)... ,(innerDia/2 + steelThick) * sin(-magnetInnerArc * pi/360 + i*pi/numberpp)); mi_selectarcsegment(2*radius - (innerDia/2 + steelThick + magnetThick) * ... cos(magnetInnerArc * pi/360 + i*pi/numberpp), (innerDia/2 + steelThick + ... magnetThick) * sin(magnetInnerArc * pi/360 + i*pi/numberpp)); mi_setarcsegmentprop(2.5, '', 0, 1); mi_selectsegment(2*radius - (innerDia/2 + steelThick + magnetThick) * ... cos(magnetInnerArc * pi/360 + i*pi/numberpp), (innerDia/2 + steelThick + ... magnetThick) * sin(magnetInnerArc * pi/360 + i*pi/numberpp)); mi_selectsegment(2*radius - (innerDia/2 + steelThick + magnetThick) * ... cos(-magnetInnerArc * pi/360 + i*pi/numberpp), (innerDia/2 + steelThick + ... magnetThick) * sin(-magnetInnerArc * pi/360 + i*pi/numberpp)); mi_setsegmentprop('', 0.5, 0, 0, 1); end mi_clearselected(); %% add block labels mi_addblocklabel(2*radius, 0); mi_addblocklabel(2*radius, outerDia*0.65); mi_clearselected(); mi_selectlabel(2*radius, 0); mi_selectlabel(2*radius, outerDia*0.65); mi_setblockprop('Air', 1, 0, '', 0, 0, 0); %% set magnetic properties mi_clearselected(); mi_addblocklabel(2*radius, outerDia/2 - steelThick - magnetThick/2); mi_addblocklabel(2*radius, innerDia/2 + steelThick + magnetThick/2); mi_selectlabel(2*radius, outerDia/2 - steelThick - magnetThick/2); mi_selectlabel(2*radius, innerDia/2 + steelThick + magnetThick/2); mi_setblockprop(magnetType, 0, 0.5, '', 90, 1, 0); mi_clearselected(); mi_selectgroup(1); mi_refreshview(); mi_copyrotate(2*radius, 0, 360/numberpp, numberpp-1); mi_selectgroup(1); mi_moverotate(2*radius, 0, 180/numberpp); mi_addblocklabel(2*radius, outerDia/2 - steelThick - magnetThick/2); mi_addblocklabel(2*radius, innerDia/2 + steelThick + magnetThick/2); mi_clearselected(); mi_selectlabel(2*radius, outerDia/2 - steelThick - magnetThick/2); mi_selectlabel(2*radius, innerDia/2 + steelThick + magnetThick/2); mi_setblockprop(magnetType, 0, 0.5, '', -90, 1, 0); mi_copyrotate(2*radius, 0, 360/numberpp, numberpp-1); mi_addblocklabel(2*radius, - outerDia/2 + steelThick/2); mi_addblocklabel(2*radius, - innerDia/2 - steelThick/2); mi_clearselected(); mi_selectlabel(2*radius, - outerDia/2 + steelThick/2); mi_selectlabel(2*radius, - innerDia/2 - steelThick/2); mi_setblockprop(steelType, 0, 0.5, '', 0, 1, 0); mi_selectgroup(1); mi_moverotate(2*radius, 0, 90/numberpp); mi_addblocklabel(2*radius, outerDia/2 - steelThick - magnetThick/2); mi_selectlabel(2*radius, outerDia/2 - steelThick - magnetThick/2); mi_setblockprop('Air', 1, 0, '', 0, 1, 0); mi_selectgroup(1); mi_moverotate(2*radius, 0, -90/numberpp); mi_refreshview(); %% draw the windings %% we set the number of bobbins at 3 x the pople pairs x reduce numberbob = 3 * numberpp * bobbinReduce ; windingInner = innerDia/2 + steelThick + magnetThick + 0.4; windingOuter = outerDia/2 - steelThick - magnetThick - 0.4; windingThick = windingOuter - windingInner; middleRad = windingInner + windingThick/2; middleCirc = middleRad * 2 * pi; coreArc = magnetInnerArc; % arc subtended by the bobbin core bobbinCore = coreArc * middleRad * pi / 180; windingArc = 180/numberbob - coreArc/2; % arc subtended by one winding bobbinMax = windingArc * windingInner * pi/180 * windingThick / wireDia^2 ; numberTurns = floor(bobbinMax * 0.78); % number of turns in the winding (empirical) mi_addnode(2*radius, windingInner); mi_addnode(2*radius, windingOuter); % first winding edge mi_addsegment(2*radius, windingInner, 2*radius, windingOuter); mi_clearselected(); mi_selectsegment(2*radius, windingInner); mi_setsegmentprop('', 0.5, 0, 0, 2); mi_clearselected(); mi_selectgroup(2); mi_moverotate(2*radius, 0, windingArc); mi_addnode(2*radius, windingInner); mi_addnode(2*radius, windingOuter); % second winding edge mi_addsegment(2*radius, windingInner, 2*radius, windingOuter); mi_clearselected(); mi_selectsegment(2*radius, windingInner); mi_setsegmentprop('', 0.5, 0, 0, 2); % outer winding arc mi_addarc(2 * radius, windingOuter, 2 * radius - windingOuter * ... sin(windingArc * pi/180),windingOuter * cos(windingArc * pi/180), ... windingArc, 0.5); % inner winding arc mi_addarc(2 * radius, windingInner, 2 * radius - windingInner * ... sin(windingArc * pi/180),windingInner * ... cos(windingArc * pi/180), windingArc, 0.5); mi_selectarcsegment(2 * radius, windingOuter); mi_selectarcsegment(2 * radius, windingInner); mi_setarcsegmentprop(1, '', 0, 2); mi_selectgroup(2); mi_mirror(2*radius, 0, 2*radius, radius); mi_selectgroup(2); mi_copyrotate(2*radius, 0, 360/numberbob, numberbob - 1); %% add block labels to the windings of phase '+A' and '-B' mi_addblocklabel(2*radius + middleRad * sin(windingArc * pi/360), middleRad); mi_addblocklabel(2*radius - middleRad * sin(windingArc * pi/360), middleRad); mi_clearselected(); mi_selectlabel(2*radius + middleRad * sin(windingArc * pi/360), middleRad); mi_setblockprop(wireType, 0, 0.5, 'A', 0, 2, numberTurns); mi_copyrotate(2*radius, 0, 360*3/numberbob, numberbob/3 -1); mi_clearselected(); mi_selectlabel(2*radius - middleRad * sin(windingArc * pi/360), middleRad); mi_setblockprop(wireType, 0, 0.5, 'B', 0, 2, -numberTurns); mi_copyrotate(2*radius, 0, 360*3/numberbob, numberbob/3 -1); %% add block labels to the windings of phase 'B' and '-C' mi_clearselected(); mi_selectgroup(2); mi_moverotate(2*radius, 0, -360/numberbob); mi_addblocklabel(2*radius + middleRad * sin(windingArc * pi/360), middleRad); mi_addblocklabel(2*radius - middleRad * sin(windingArc * pi/360), middleRad); mi_clearselected(); mi_selectlabel(2*radius + middleRad * sin(windingArc * pi/360), middleRad); mi_setblockprop(wireType, 0, 0.5, 'B', 0, 2, numberTurns); mi_copyrotate(2*radius, 0, 360*3/numberbob, numberbob/3 -1); mi_clearselected(); mi_selectlabel(2*radius - middleRad * sin(windingArc * pi/360), middleRad); mi_setblockprop(wireType, 0, 0.5, 'C', 0, 2, -numberTurns); mi_copyrotate(2*radius, 0, 360*3/numberbob, numberbob/3 -1); %% add block labels to the windings of phase 'C' and '-A' mi_clearselected(); mi_selectgroup(2); mi_moverotate(2*radius, 0, -360/numberbob); mi_addblocklabel(2*radius + middleRad * sin(windingArc * pi/360), middleRad); mi_addblocklabel(2*radius - middleRad * sin(windingArc * pi/360), middleRad); mi_clearselected(); mi_selectlabel(2*radius + middleRad * sin(windingArc * pi/360), middleRad); mi_setblockprop(wireType, 0, 0.5, 'C', 0, 2, numberTurns); mi_copyrotate(2*radius, 0, 360*3/numberbob, numberbob/3 -1); mi_clearselected(); mi_selectlabel(2*radius - middleRad * sin(windingArc * pi/360), middleRad); mi_setblockprop(wireType, 0, 0.5, 'A', 0, 2, -numberTurns); mi_copyrotate(2*radius, 0, 360*3/numberbob, numberbob/3 -1); mi_selectgroup(2); mi_moverotate(2*radius, 0, -180/numberbob); mi_refreshview(); %%---------------------------------------------------------------------------------------------- mi_saveas('Windgen-test.FEM'); outFile = 'FEMMresult.txt'; fid = fopen(outFile,'w'); % if outFile==nil print('Failed to create result file: ', why) endif %% to remove uncertainty, write out a header in the results file fprintf(fid, 'Coreless wind turbine generator by Jose B. Almeida, May 2015. \n'); fprintf(fid, 'OctaveFEMM file: Arc_magnets.m \n'); fprintf(fid, 'FEM file: Windgen-test.FEM \n \n'); fprintf(fid, '# Number of pole pairs = %d \n', numberpp); fprintf(fid, '# Magnet type = %s \n', magnetType); fprintf(fid, '# Outer diameter = %u %s \n', outerDia, units); fprintf(fid, '# Inner diameter = %u %s \n', innerDia, units); fprintf(fid, '# Magnet l x t = %u x %u %s \n',magnetLength, magnetThick, units); fprintf(fid, '# Outer magnet arc = %f degrees \n', magnetOuterArc); fprintf(fid, '# Outer magnet radii = %f, %f %s \n', outerDia/2 - steelThick, outerDia/2 - steelThick - magnetThick, units); fprintf(fid, '# Inner magnet arc = %f degrees \n', magnetInnerArc); fprintf(fid, '# Inner magnet radii = %f, %f %s \n', innerDia/2 + steelThick + magnetThick, innerDia/2 + steelThick, units); fprintf(fid, '# Number of bobbins = %u \n', numberbob); fprintf(fid, '# Bobbin thickness = %f %s \n', windingThick, units); fprintf(fid, '# Width of bobbin core = %u %s \n', bobbinCore, units); fprintf(fid, '# Wire diameter = %f %s \n', wireDia, units); fprintf(fid, '# Number of turns = %u \n', numberTurns); %%---------------------------------------------------------------------------------------------- % start processing mi_analyze; mi_loadsolution; % opens the solution window %% pole pitch span contour for interactive analysis %% set the window to a nice size for the problem mo_zoom(radius*1.2, 0.0, radius*1.5, radius*0.75); mo_showdensityplot(1,0,2.0,0.0,'mag'); %%Shows the flux density mo_seteditmode('contour'); mo_addcontour(2*radius, middleRad); mo_addcontour(2*radius + middleRad * sin(4 * pi/numberpp), middleRad * ... cos(4 * pi/numberpp)); mo_bendcontour(-720/numberpp, 1); %% evaluate normal induction for plotting and saving n = 100; for k = 1:n alpha(k) = k/n * 4 * pi/numberpp; point = mo_getb(2*radius + middleRad * sin(alpha(k)), middleRad * ... cos(alpha(k))); induction(k) = point(1)*sin(alpha(k)) + point(2)* cos(alpha(k)); end %% Sinusoid for comparison maxInduction = max(induction); n = 100; for k = 1:n sinus(k) = maxInduction*sin(alpha(k)*numberpp); end alphaDegree = alpha * 180/pi; plot(alphaDegree,induction,alphaDegree,sinus); title('B radial (blue) and Sinusoid (red) comparison','fontsize',12); xlabel('Rotor Angle, Degrees'); ylabel('B-radial, Tesla'); %axis([0, alphaDegree(n)]); fprintf(fid, '\n \n# Induction amplitude = %f Tesla \n', maxInduction); fclose(fid);