Verwendung von poc-Files in individuellen Berechnungsabläufen¶
Skript-Datei
--------------------------------------------------------------------------------
-- General Settings ------------------------------------------------------------
--------------------------------------------------------------------------------
exit_on_error = false -- Action after error
exit_on_end = true -- Action after script termination
verbosity = 3 -- Level of feedback
--------------------------------------------------------------------------------
-- Model Definitions -----------------------------------------------------------
--------------------------------------------------------------------------------
Q = 12 -- Number of slots
P = 10 -- Number of poles
Da = 100 -- Stator outer diameter
Di = 55 -- Bore diameter
s = 5 -- Slot width
ag = 1 -- Air-gap width
bz = 7 -- Tooth width
h1 = 1.5 -- Tooth tip height 1
h2 = 2 -- Tooth tip height 2
hj = 8 -- Yoke height
hrs = 6 -- Rotor back iron height
hm = 3 -- Magnet height
bm = 11.205 -- Magnet width
alpham = 0.88 -- Polbedeckung
ls = 150 -- Stack length
mtype = 2 -- Type of magnet: 1 = ring segment, 2 = cuboid
Qm = Q/2 -- Number of slots in model
Pm = P/2 -- Number of poles in model
urs = 1000 -- Permeability stator iron
urr = 1000 -- Permeability rotor back iron
Br = 1.2 -- Permanent magnet remanence
urm = 1.05 -- Permanent magnet permeability
Nc = 10 -- Number of turns (coil)
--------------------------------------------------------------------------------
-- Model Generation ------------------------------------------------------------
--------------------------------------------------------------------------------
new_model_force("example","PMSM")
global_unit("mm") -- Global unit (m, cm, mm)
pickdist(0.001) -- Snap distance
blow_up_wind(0,0,55,55) -- Window size
cosys("cartes")
------------
-- Stator --
------------
-- Calculation of characteristic points
x,y = {}, {}
for i=1, 15 do
x[i]=0
y[i]=0
end
x[1],y[1] = pd2c(Da/2,0)
x[2],y[2] = pd2c(Da/2,180/Q)
x[3] = Di/2*math.cos(math.asin(s/Di))
y[3] = s/2
x[4] = x[3]+h1
y[4] = s/2
x[5],y[5] = pd2c(Di/2,180/Q)
x[6] = Di/2+h1+h2;
y[6] = y[5]/x[5]*x[6]-bz/2/math.cos(pi/Di)
x[7] = Da/2-hj;
y[7] = y[5]/x[5]*x[7]-bz/2/math.cos(pi/Di)
x[8] = x[7]
x[9],y[9] = pd2c(vlen(x[4],y[4]),180/Q)
x[10] = (y[6]+x[5]/y[5]*x[6])/(y[5]/x[5]+x[5]/y[5])
y[10] = y[5]/x[5]*x[10]
x[11] = (y[7]+x[5]/y[5]*x[7])/(y[5]/x[5]+x[5]/y[5])
y[11] = y[5]/x[5]*x[11]
x[12] = Di/2
x[13] = x[4]
x[14] = Di/2-ag/3
x[15],y[15] = pd2c(Di/2-ag/3,180/Q)
-- Node chain generation
agnp = 1 -- Node pitch (degree) in cirumferential direction
ndt(ag) -- reference node distance
nc_circle(x[14],y[14],x[15],y[15],360/Q/2/agnp+1)
nc_circle(x[1],y[1],x[2],y[2],0)
nc_circle(x[13],y[13],x[4],y[4],0)
nc_circle(x[3],y[3],x[5],y[5],0)
nc_line(x[3],y[3],x[4],y[4],0)
nc_line_cont(x[6],y[6],0)
nc_line_cont(x[7],y[7],0)
nc_line_cont(x[8],y[8],0)
nc_line(x[12],y[12],x[13],y[13],0)
nc_line_cont(x[8],y[8],0)
nc_line_cont(x[1],y[1],0)
nc_line(x[14],y[14],x[12],y[12],0)
nc_line(x[15],y[15],x[5],y[5],0)
nc_line_cont(x[9],y[9],0)
nc_line_cont(x[10],y[10],0)
nc_line_cont(x[11],y[11],0)
nc_line_cont(x[2],y[2],0)
-- Meshing
mesh.con1 = 0.1 -- Mesh control (growth of element size)
create_mesh_se(Da/2-hj/2,0+hj/2)
create_mesh_se((Da+Di)/4,s/4)
create_mesh_se(Di/2+h1/2,s/4)
-- Definition of subregions
def_new_subreg(Da/2-hj/2,0+hj/2,"Stator",11)
-- Mirror, rotate and copy regions
mirror_nodechains(x[2],y[2],x[15],y[15])
x0,y0 = pd2c(Di/2-ag/3,0)
x1,y1 = pd2c(Da/2,0)
x2,y2 = pd2c(Da/2,360/Q)
x3,y3 = pd2c(Di/2-ag/3,360/Q)
rotate_copy_nodechains(x0,y0,x1,y1,x2,y2,x3,y3,Qm-1)
-----------
-- Rotor --
-----------
-- Calculation of characteristic points
x[1],y[1] = pd2c(Di/2-ag*2/3,0)
x[2],y[2] = pd2c(Di/2-ag*2/3,360/P)
x[3],y[3] = pd2c(Di/2-ag,0)
x[4],y[4] = pd2c(Di/2-ag,360/P)
x[5],y[5] = pd2c(Di/2-ag-hm,0)
x[6],y[6] = pd2c(Di/2-ag-hm,360/P)
x[7],y[7] = pd2c(Di/2-ag-hm-hrs,0)
x[8],y[8] = pd2c(Di/2-ag-hm-hrs,360/P)
x[8],y[8] = pd2c(Di/2-ag-hm-hrs,360/P)
if (mtype==1) then -- ring segment
x[9],y[9] = pd2c(Di/2-ag,360/P*(1-alpham)*0.5)
x[10],y[10] = pd2c(Di/2-ag,360/P*(1+alpham)*0.5)
x[11],y[11] = pd2c(Di/2-ag-hm,360/P*(1-alpham)*0.5)
x[12],y[12] = pd2c(Di/2-ag-hm,360/P*(1+alpham)*0.5)
else -- cuboid
R = Di/2-ag-hm
alpha = math.asin(bm/(2*R))
x[11],y[11] = pd2c(R,360/P*0.5-alpha/math.pi*180)
x[12],y[12] = pd2c(R,360/P*0.5+alpha/math.pi*180)
R = Di/2-ag
alpha = math.asin(bm/(2*R))
x[9],y[9] = pd2c(R,360/P*0.5-alpha/math.pi*180)
x[10],y[10] = pd2c(R,360/P*0.5+alpha/math.pi*180)
end
-- Node chain generation
ndt(0.75) -- reference node distance
nc_circle(x[1],y[1],x[2],y[2],360/P/agnp+1)
nc_circle(x[3],y[3],x[9],y[9],0)
nc_circle(x[9],y[9],x[10],y[10],0)
nc_circle(x[10],y[10],x[4],y[4],0)
nc_circle(x[5],y[5],x[11],y[11],0)
if (mtype==1) then
nc_circle(x[11],y[11],x[12],y[12],0)
else
nc_line(x[11],y[11],x[12],y[12],0)
nc_line(x[9],y[9],x[10],y[10],0)
end
nc_circle(x[12],y[12],x[6],y[6],0)
nc_circle(x[7],y[7],x[8],y[8],0)
nc_line(x[7],y[7],x[5],y[5],0)
nc_line_cont(x[3],y[3],0)
nc_line_cont(x[1],y[1],0)
nc_line(x[8],y[8],x[6],y[6],0)
nc_line_cont(x[4],y[4],0)
nc_line_cont(x[2],y[2],0)
nc_line(x[11],y[11],x[9],y[9],0)
nc_line(x[12],y[12],x[10],y[10],0)
-------------
-- Meshing --
-------------
create_mesh_se(pd2c(Di/2-ag*5/6,180/P))
create_mesh_se(pd2c(Di/2-ag-hm/2,180/P))
create_mesh_se(pd2c(Di/2-ag-hm-hrs/2,180/P))
create_mesh_se(pd2c(Di/2-ag-hm/2,360/P*(1-alpham)*0.25))
create_mesh_se(pd2c(Di/2-ag-hm/2,360/P*(3+alpham)*0.25))
if (mtype == 2) then
R = Di/2-ag-hm
alpha = math.asin(bm/(2*R))
Rx = R*math.cos(alpha)
create_mesh_se(pd2c(Di/2-ag-(R-Rx)/2,360/P*0.5))
end
-- Definition of subregions
def_new_subreg(Di/2-ag-hm-hrs/2,ag,"Rückschluss",11)
-- Mirror, rotate and copy regions
rotate_copy_nodechains(x[7],y[7],x[1],y[1],x[2],y[2],x[8],y[8],Pm-1)
-- Air gap
x0,y0 = pd2c(Di/2-ag*2/3,0)
x1,y1 = pd2c(Di/2-ag/3,0)
nc_line(x0,y0,x1,y1,0)
x0,y0 = pd2c(Di/2-ag*2/3,360*Pm/P)
x1,y1 = pd2c(Di/2-ag/3,360*Pm/P)
nc_line(x0,y0,x1,y1,0)
create_mesh_se(Di/2-ag/2,ag)
for i=1, 12 do -- output points and point numbers for debugging purposes
point(x[i],y[i],"red",".")
text(x[i],y[i]+1,tostring(i),"red",0.3)
end
--debug()
-------------------------
-- Boundary Conditions --
-------------------------
x0,y0 = pd2c(Di/2-ag-hm-hj,0)
x1,y1 = pd2c(Da/2,0)
x2,y2 = pd2c(Di/2-ag-hm-hj,360*Pm/P)
x3,y3 = pd2c(Da/2,360*Pm/P)
def_bcond_vpo(x1,y1,x3,y3)
def_bcond_vpo(x2,y2,x0,y0)
def_bcond(x3,y3,x2,y2,x0,y0,x1,y1,"neg")
--------------
-- Windings --
--------------
tauq = 360/Q -- Nutteilungswinkel
Rq = (Di/2+Da/2-hj)/2 -- mittlerer Nutradius
x,y = pd2c(Rq,tauq/4)
wkey = def_new_wdg(x,y,"cyan","Strang 1",Nc,0.0,"wo")
x,y = pd2c(Rq,tauq-tauq/4)
add_to_wdg(x,y,"wsamekey","wi","wser")
x,y = pd2c(Rq,tauq+tauq/4)
add_to_wdg(x,y,"wsamekey","wi","wser")
x,y = pd2c(Rq,2*tauq-tauq/4)
add_to_wdg(x,y,"wsamekey","wo","wser")
x,y = pd2c(Rq,2*tauq+tauq/4)
def_new_wdg(x,y,"magenta","Strang 2",Nc,0.0,"wi")
x,y = pd2c(Rq,3*tauq-tauq/4)
add_to_wdg(x,y,"wsamekey","wo","wser")
x,y = pd2c(Rq,3*tauq+tauq/4)
add_to_wdg(x,y,"wsamekey","wo","wser")
x,y = pd2c(Rq,4*tauq-tauq/4)
add_to_wdg(x,y,"wsamekey","wi","wser")
x,y = pd2c(Rq,4*tauq+tauq/4)
def_new_wdg(x,y,"yellow","Strang 3",Nc,0.0,"wo")
x,y = pd2c(Rq,5*tauq-tauq/4)
add_to_wdg(x,y,"wsamekey","wi","wser")
x,y = pd2c(Rq,5*tauq+tauq/4)
add_to_wdg(x,y,"wsamekey","wi","wser")
x,y = pd2c(Rq,6*tauq-tauq/4)
add_to_wdg(x,y,"wsamekey","wo","wser")
---------------
-- Materials --
---------------
-- Stator und rotor back iron
def_mat_fm(Da/2-hj/2,ag,urs,100)
def_mat_fm(Di/2-ag-hm-hrs/2,ag,urr,100)
-- Permanent magnets
if (mtype == 2) then
for i=0, Pm/2 do
alpha = 360/P*(2*i+1)-180/P
x,y = pd2c(Di/2-ag-hm/2,alpha)
def_mat_pm(x,y,"red",Br,urm,alpha,"parallel",100)
end
for i=1, Pm/2 do
alpha = 360/P*2*i-180/P
x,y = pd2c(Di/2-ag-hm/2,alpha)
def_mat_pm(x,y,"green",Br,urm,alpha+180,"parallel",100)
end
else
for i=0, Pm/2 do
x,y = pd2c(Di/2-ag-hm/2,360/P*(2*i+1)-180/P)
def_mat_pm(x,y,"red",Br,urm,0,"radial",100)
end
for i=1, Pm/2 do
x,y = pd2c(Di/2-ag-hm/2,360/P*2*i-180/P)
def_mat_pm(x,y,"green",Br,urm,180,"radial",100)
end
end
---------------------------------------
-- Set basic machine and model data ---
---------------------------------------
m.tot_num_slot = Q -- number of Slots
m.num_slots = Qm -- number of Slots simulated
m.num_poles = P -- number of Poles 2p (>= 2)
m.npols_gen = Pm -- number of Poles simulated (>= 1)
m.arm_length = ls -- stack length
pre_models("basic_modpar");
------------------------------
-- Generator poc file for ---
------------------------------
cosys("polar") -- refeference system has to be polar (r-phi)
pre_models("gen_pocfile");
--------------------------------------------------------------------------------
-- Calculations ----------------------------------------------------------------
--------------------------------------------------------------------------------
-- define some general data related to following multi calculations
Rag = 0.5*(Di-ag) -- average air gap radius
phi = 0.0 -- initial position
dphi = 1.0 -- angular displacement incement
NRot = 720.0/P/dphi+1 -- number of movement steps (including first and last)
current = 30 -- current (peak)
angl_i_up = -20 -- current angle
verbosity = 1
outputfile=io.open("output.txt","w+")
outputfile:write("# phi [°] Psi_wk1 [Vs] T [Nm]\n");
px, pyT, pyP = {}, {}, {} -- define arrays for data to plot later
-- read and process the poc file
read_poc_file("example_10p.poc",NRot,current,angl_i_up)
-- initialize the movement
rotate({
airgap = Rag, -- air gap radius
region = "inside", -- region to rotate
mode = "save", -- save initial model state
quant = dphi -- quantization of rotation angle
})
phi = 0.0; -- set instantanious position variable to zero
wkey = 1 -- first winding
for i=1,NRot do -- loop over all positions
phi, dphi = update_curr_pos(i) -- update winding currents for position i
-- this function returns the current position
-- and the position increment (used later to move next)
calc_field_single({ -- single field calculation
maxit = 1, -- max number of iterations
maxcop = 0.005, -- max change of permeability
permode = "act" -- permeability mode (use instantaneous permeability)
})
Psi = flux_winding_wk(wkey)*ls*P/Pm -- flux linkage evaluation in winding
m.coord_x1, m.coord_y1 = pd2c(Rag,0.0) -- torque calculation
m.coord_x2, m.coord_y2 = pd2c(Rag,359.5) -- (not very comfortable yet, will be changed in future)
post_models("force_torque","FT")
T = FT[3]*P/Pm -- calculated overall torque based on modelled part
-- output results to console and to ascii file
printf(" > phi = %g deg, Psi = %g Vs, M = %g Nm",phi,Psi,T)
outputfile:write(string.format("%7.3f %9.6f %9.6f\n",phi,Psi,T));
if (phi<0.5*360.0*Pm/P) then -- process movement steps
rotate({
angle = phi+dphi, -- angle of rotation
mode = "abs", -- type of movement (here: absolute angle)
check = "auto", -- type of internal checks to be performed
})
else
rotate({ -- avoid rotate outside the model boundaries
angle = -360.0*Pm/P+phi+dphi, -- shift rotor to equivalent angular position inside
})
end
px[i] = phi -- save instantaneous position,
pyT[i] = T -- torque
pyP[i] = Psi -- and flux linkage for plotting
if (i>2) then -- display diagrams after the second movement steps
plot({
x = px, -- torque diagram
y = pyT,
pos = 21,
xmin = 0,
xmax = 360/P*2,
ymin = 0,
xlabel = "Rotor position [Dgr]",
ylabel = "Torque [Nm]",
color = "red"
})
plot({
x = px, -- flux linkage diagram
y = pyP,
pos = 22,
xmin = 0,
xmax = 360/P*2,
xlabel = "Rotor position [Dgr]",
ylabel = "Flux linkage [Vs]",
color = "blue"
})
end
end
io.close(outputfile)
save_metafile("example_pocfile.eps") -- save graphics output to file
rotate({
mode = "reset" -- reset model to initial state and discard changes
})
save_model("close")
