Расчёт по полному угловому периоду и спектральный анализ

"""Встроенный Python FEMM5. Три учебные геометрии на основе Katona, независимый адаптер.
Новые серии: 32 положения ХХ и 32 положения нагрузки на зубцовый период.
Все токи пиковые. Nф=40 последовательных витков, Iq=5А, d=40мм.
"""
import math,json,os,sys
assert sys.implementation.name=='micropython'

# FEMM CLI и команда запуска сохранённого файла передают его путь в sys.argv[0].
# Локальные модели/таблицы ищутся возле программы, а не в случайной рабочей папке.
import os
_saved_file=globals().get('__file__') or (sys.argv[0] if sys.argv else '')
if not _saved_file or not os.path.isfile(_saved_file):
 raise RuntimeError('Запустите сохранённый файл из полного распакованного комплекта командой FEMM Run Script; не вставляйте его частями в консоль.')
os.chdir(os.path.dirname(os.path.abspath(_saved_file)))

try:os.mkdir('katona-output')
except OSError:pass
def close():
 for n in ['mo_close','mi_close']:
  try:globals()[n]()
  except Exception:pass
def labels(path):
 with open(path,'r') as f:lines=f.read().split('\n')
 names=[];rows=[];inlabels=False;count=0
 for line in lines:
  if '<BlockName>' in line:names.append(line.split('"')[1])
  if '[NumBlockLabels]' in line:inlabels=True;count=int(line.split('=')[1]);continue
  if inlabels and len(rows)<count:rows.append(line.split())
 return names,rows
def load(case):
 path='katona-'+str(case)+'-course.fem';close();opendocument(path);mi_probdef(0,'millimeters','planar',1e-8,40,25)
 materials,rows=labels(path);counts=[0,0,0]
 for row in rows:
  name=materials[int(row[2])-1]
  if name in ['A','B','C','a','b','c']:counts['ABC'.index(name.upper())]+=1
 assert len(set(counts))==1
 # Каждая сторона катушки содержит N, число последовательных витков фазы = сумма |N|/2.
 turns=80/counts[0]
 for name in ['A','B','C']:mi_addcircprop('Фаза '+name,0,1)
 mi_addmaterial('Идеальная обмотка',1,1,0,0,0)
 for row in rows:
  mat=materials[int(row[2])-1]
  if mat in ['A','B','C','a','b','c']:
   x,y=float(row[0]),float(row[1]);mi_selectlabel(x,y);mi_setblockprop('Идеальная обмотка',0,1,'Фаза '+mat.upper(),0,0,turns if mat.isupper() else -turns);mi_clearselected()
 # Группа ротора заранее назначена по радиусам обоих концов каждой линии/дуги.
 # Выбор окружностью опасен для длинной внешней дуги, чья хорда проходит через ротор.
 # Воздушная метка зазора не входит в интегрируемое тело ротора.
 for row in rows:
  mat=materials[int(row[2])-1]
  if mat=='Воздух':
   x,y=float(row[0]),float(row[1]);mi_selectlabel(x,y);mi_setblockprop(mat,1,0,'',0,0,0);mi_clearselected()
 return {'coil_side_count_per_phase':counts,'turns_per_side':turns,'series_turns_per_phase':40}
def sample(case,name,theta,beta,offset):
 load(case);p=2 if case<3 else 4
 mi_selectgroup(2);mi_moverotate(0,0,theta);mi_clearselected()
 angle=math.radians(p*theta)+offset
 currents=[0.,0.,0.] if name=='open' else [-5*math.sin(angle-k*2*math.pi/3) for k in range(3)]
 for k,n in enumerate(['A','B','C']):mi_modifycircprop('Фаза '+n,1,currents[k])
 path='katona-output/case-'+str(case)+'-'+name+'-'+str(round(theta,4))+'.fem';mi_saveas(path);mi_analyze(1);mi_loadsolution()
 psi=[float(mo_getcircuitproperties('Фаза '+n)[2]) for n in ['A','B','C']]
 assert all(math.isfinite(v) for v in psi),'Нефинитное потокосцепление: проверьте тип проводника и проводимость'
 mo_clearblock();mo_groupselectblock(2);T=float(mo_blockintegral(22));mo_clearblock();row={'mechanical_angle_deg':theta,'currents_A':currents,'flux_Wb':psi,'T_Nm':T,'elements':mo_numelements()};close();return row
N=2 if len(sys.argv)>1 and sys.argv[1]=='quick' else 32
out={'source':'Katona three documented FEM geometries; independent FEMM5 embedded-Python adapter, not the original external py2femm program.','scope':'Сравнение топологий с одинаковыми Nф=40, Iq=5А и глубиной40мм. Пазовые области идеализированы как заданные токовые обмотки с нулевой проводимостью; потери меди не решаются. Геометрии исходника, число полюсов, площадь меди и длина лобовых частей различаются: это не изолированный опыт укорочения шага и не сравнение КПД. Расчётные числа не являются числами публикации Katona.','angles_per_period':N,'cases':[]}
for case in [1,2,3]:
 meta=load(case);close();zero=sample(case,'open',0,0,0)
 alpha=2/3*(zero['flux_Wb'][0]-.5*zero['flux_Wb'][1]-.5*zero['flux_Wb'][2]);beta=2/3*math.sqrt(3)/2*(zero['flux_Wb'][1]-zero['flux_Wb'][2]);offset=math.atan2(beta,alpha)
 period=30 if case<3 else 15;series={'open':[],'load':[]}
 for state in ['open','load']:
  for k in range(N):series[state].append(sample(case,state,k*period/N,0,offset))
 stats={}
 for state,rows in series.items():
  v=[r['T_Nm'] for r in rows];coarse=v[::2];stats[state]={'mean_T_Nm':sum(v)/len(v),'ripple_pp_Nm':max(v)-min(v),'coarse_16_mean_Nm':sum(coarse)/len(coarse),'coarse_16_ripple_pp_Nm':max(coarse)-min(coarse)}
  if state=='open':
   vd=vq=0
   for row in rows:
    angle=math.radians((2 if case<3 else 4)*row['mechanical_angle_deg'])+offset
    vd+=2/3*sum(row['flux_Wb'][k]*math.cos(angle-k*2*math.pi/3) for k in range(3))
    vq-=2/3*sum(row['flux_Wb'][k]*math.sin(angle-k*2*math.pi/3) for k in range(3))
   stats[state]['pm_linkage_fundamental_projection_Wb']=math.sqrt(vd*vd+vq*vq)/len(rows)
 out['cases'].append({'case':case,'p':2 if case<3 else 4,'mechanical_period_deg':period,'electrical_offset_rad':offset,'winding':meta,'series':series,'statistics':stats})
with open('katona-control-result.json' if N<32 else 'katona-result.json','w') as f:json.dump(out,f)
print('Katona comparison',[(c['case'],c['statistics']) for c in out['cases']])

"""MicroPython FEMM5. Многослойная 2D-оценка по фактической кривой ХХ.
Срезы независимы; это постобработка, а не трёхмерное полевое решение.
"""
import math,json,sys
assert sys.implementation.name=='micropython'

# FEMM CLI и команда запуска сохранённого файла передают его путь в sys.argv[0].
# Локальные модели/таблицы ищутся возле программы, а не в случайной рабочей папке.
import os
_saved_file=globals().get('__file__') or (sys.argv[0] if sys.argv else '')
if not _saved_file or not os.path.isfile(_saved_file):
 raise RuntimeError('Запустите сохранённый файл из полного распакованного комплекта командой FEMM Run Script; не вставляйте его частями в консоль.')
os.chdir(os.path.dirname(os.path.abspath(_saved_file)))

with open('spm-open-compact.json') as f:data=json.load(f)
rows=data['rows'];N=len(rows);step=data['step_mechanical_deg'];values=[r[1] for r in rows]
shift=round(7.5/step);assert abs(shift*step-7.5)<1e-10
skew=[.5*(values[(k-shift)%N]+values[(k+shift)%N]) for k in range(N)]
def rms(v):
 mean=sum(v)/len(v);return math.sqrt(sum((x-mean)**2 for x in v)/len(v))
result={'source':data['source'],'scope':'Линейное объединение независимых двумерных срезов по исходной полевой угловой кривой; новые поля и осевые потоки не рассчитываются. Исходная кривая ХХ, поэтому о снижении среднего нагрузочного момента по ней судить нельзя.','N':N,'offsets_mechanical_deg':[-7.5,7.5],'weights':[.5,.5],'base_cogging_pp_Nm':max(values)-min(values),'skew_cogging_pp_Nm':max(skew)-min(skew),'base_ripple_rms_Nm':rms(values),'skew_ripple_rms_Nm':rms(skew),'ideal_fundamental_EMF_ratio_two_slices':math.cos(2*math.radians(7.5)),'ideal_continuous_30deg_EMF_ratio':math.sin(math.pi/6)/(math.pi/6),'rows':[[rows[k][0],values[k],skew[k]] for k in range(N)]}
with open('skew-result.json','w') as f:json.dump(result,f)
print('Skew pp',result['base_cogging_pp_Nm'],result['skew_cogging_pp_Nm'])

"""FEMM5 MicroPython. Независимый контролируемый опыт изменения дуги PM.
Геометрия строится заново; оригинальные файлы курса не изменяются.
"""
import math,json,os,sys
assert sys.implementation.name=='micropython'

# FEMM CLI и команда запуска сохранённого файла передают его путь в sys.argv[0].
# Локальные модели/таблицы ищутся возле программы, а не в случайной рабочей папке.
import os
_saved_file=globals().get('__file__') or (sys.argv[0] if sys.argv else '')
if not _saved_file or not os.path.isfile(_saved_file):
 raise RuntimeError('Запустите сохранённый файл из полного распакованного комплекта командой FEMM Run Script; не вставляйте его частями в консоль.')
os.chdir(os.path.dirname(os.path.abspath(_saved_file)))

MU0=4e-7*math.pi
nodes={};edges=set();arcs=set();global_angles=[];pm_label_specs=[]
try:os.mkdir('magnet-arc-output')
except OSError:pass
def node(x,y):
    key=(round(x,7),round(y,7))
    if key not in nodes:mi_addnode(x,y);nodes[key]=1

def line(x1,y1,x2,y2):
    node(x1,y1);node(x2,y2)
    key=tuple(sorted(((round(x1,7),round(y1,7)),(round(x2,7),round(y2,7)))))
    if key not in edges:mi_addsegment(x1,y1,x2,y2);edges.add(key)

def rect(x1,y1,x2,y2):
    line(x1,y1,x2,y1);line(x2,y1,x2,y2);line(x2,y2,x1,y2);line(x1,y2,x1,y1)

def polar(r,a):
    return r*math.cos(math.radians(a)), r*math.sin(math.radians(a))

def arc(r,a,b):
    a=a%360; span=(b-a)%360
    if span<1e-7:return
    cuts=[a]+sorted([a+(t-a)%360 for t in global_angles if 1e-6<(t-a)%360<span-1e-6])+[a+span]
    for left,right in zip(cuts[:-1],cuts[1:]):
        x1,y1=polar(r,left);x2,y2=polar(r,right);node(x1,y1);node(x2,y2)
        key=(round(r,7),round(left%360,7),round(right%360,7))
        if key not in arcs:mi_addarc(x1,y1,x2,y2,right-left,4);arcs.add(key)

def circle(r):
    angles=sorted(set([0,90,180,270]+[round(a%360,7) for a in global_angles]))
    for a,b in zip(angles,angles[1:]+[angles[0]+360]):arc(r,a,b)

def sector(r1,r2,a,b):
    arc(r1,a,b);arc(r2,a,b)
    x1,y1=polar(r1,a);x2,y2=polar(r2,a);line(x1,y1,x2,y2)
    x1,y1=polar(r1,b);x2,y2=polar(r2,b);line(x1,y1,x2,y2)

def label(x,y,material,circuit='',turns=0,mesh=2,group=0,magdir=0):
    mi_addblocklabel(x,y);mi_selectlabel(x,y)
    mi_setblockprop(material,0,mesh,circuit,magdir,group,turns);mi_clearselected()
    if material=='Magnet':
        pm_label_specs.append({'x':x,'y':y,'direction':magdir,
                               'numeric_angle_deg':0 if isinstance(magdir,str) else magdir})

def boundary_segment(x,y):
    mi_selectsegment(x,y);mi_setsegmentprop('A0',0,1,0,0);mi_clearselected()

def setup(frequency=0,kind='planar',depth=40):
    global nodes,edges,arcs,global_angles,pm_label_specs
    nodes={};edges=set();arcs=set();global_angles=[];pm_label_specs=[]
    newdocument(0);mi_probdef(frequency,'millimeters',kind,1e-8,depth,25)
    mi_addmaterial('Air',1,1,0,0,0)
    mi_addmaterial('Winding',1,1,0,0,0)
    mi_addmaterial('LinearSteel',1000,1000,0,0,0)
    mi_addmaterial('NonlinearSteel',1,1,0,0,0)
    for H,B in [(0,0),(100,.25),(200,.5),(400,.9),(800,1.2),(2000,1.45),(5000,1.6),(10000,1.7)]:mi_addbhpoint('NonlinearSteel',B,H)
    mi_addmaterial('Magnet',1.05,1.05,1.2/(MU0*1.05),0,0)
    mi_addboundprop('A0',0,0,0,0,0,0,0,0,0)

def outer_box(size):
    rect(-size,-size,size,size)
    for x,y in [(0,-size),(size,0),(0,size),(-size,0)]:boundary_segment(x,y)
    label(size*.85,size*.85,'Air',mesh=size/20)

def stator(p=2,current=2,angle=0):
    q=12//(6*p);mapping=[('A',1),('C',-1),('B',1),('A',-1),('C',1),('B',-1)]
    for c,a in [('A',0),('B',-120),('C',120)]:mi_addcircprop(c,current*math.cos(math.radians(angle+a)),1)
    circle(60);circle(100);label(0,94,'LinearSteel',mesh=3,group=1)
    for k in range(12):
        center=k*30;sector(60,85,center-7,center+7);x,y=polar(73,center);c,sign=mapping[(k//q)%6];label(x,y,'Winding',c,sign*40,mesh=2,group=1)
def close():
 for n in ['mo_close','mi_close']:
  try:globals()[n]()
  except Exception:pass
def one(half,theta,current,gamma,name):
 close();setup(depth=80)
 global global_angles
 global_angles=[0,90,180,270]
 for k in range(12):global_angles.extend([k*30-7,k*30+7])
 for k in range(4):global_angles.extend([theta+k*90-half,theta+k*90+half])
 outer_box(170);stator(2,current,gamma);circle(18);circle(48);circle(57)
 label(*polar(32,theta+45),'LinearSteel',mesh=2,group=2)
 for k in range(4):
  center=theta+k*90;sector(48,57,center-half,center+half)
  label(*polar(52,center),'Magnet',mesh=1.5,group=2,magdir='theta' if k%2==0 else 'theta+180')
  label(*polar(52,center+45),'Air',mesh=1)
 label(0,0,'Air',mesh=2);label(*polar(59,17),'Air',mesh=.35)
 mi_selectcircle(0,0,58.01,4);mi_setgroup(2);mi_clearselected()
 mi_saveas('magnet-arc-output/'+name+'.fem');mi_analyze(1);mi_loadsolution()
 psi=[float(mo_getcircuitproperties(n)[2]) for n in ['A','B','C']]
 assert all(math.isfinite(v) for v in psi)
 mo_clearblock();mo_groupselectblock(2);T=float(mo_blockintegral(22));mo_clearblock()
 return {'theta_mechanical_deg':theta,'flux_Wb':psi,'T_Nm':T,'elements':mo_numelements()}
N=4 if len(sys.argv)>1 and sys.argv[1]=='quick' else 32
results=[]
for half in [40.,36.]:
 first=one(half,0,0,0,'cal-'+str(half));psi=first['flux_Wb'];alpha=2/3*(psi[0]-.5*psi[1]-.5*psi[2]);beta=math.sqrt(3)/3*(psi[1]-psi[2]);offset=math.atan2(beta,alpha)
 series={};stats={}
 for state in ['open','load']:
  rows=[]
  for k in range(N):
   theta=30*k/N;gamma=2*theta+math.degrees(offset)+90
   rows.append(one(half,theta,0 if state=='open' else 5,gamma,('control-' if N<32 else '')+str(half)+'-'+state+'-'+str(k)))
  series[state]=rows;v=[r['T_Nm'] for r in rows];stats[state]={'mean_T_Nm':sum(v)/len(v),'ripple_pp_Nm':max(v)-min(v),'coarse16_mean_Nm':sum(v[::2])/(N//2),'coarse16_ripple_pp_Nm':max(v[::2])-min(v[::2])}
  if state=='open':
   pd=pq=0.
   for row in rows:
    angle=math.radians(2*row['theta_mechanical_deg'])+offset
    pd+=2/3*sum(row['flux_Wb'][j]*math.cos(angle-j*2*math.pi/3) for j in range(3))
    pq-=2/3*sum(row['flux_Wb'][j]*math.sin(angle-j*2*math.pi/3) for j in range(3))
   amp=math.sqrt(pd*pd+pq*pq)/N;stats[state]['pm_flux_fundamental_Wb']=amp;stats[state]['phase_emf_rms_V_at_1000rpm']=2*2*math.pi*1000/60*amp/math.sqrt(2)
 results.append({'full_magnet_arc_deg':2*half,'magnet_volume_relative':half/40,'series':series,'statistics':stats})
close()
out={'scope':'Собственная SPM Q12/2p4, p2, глубина80мм, PM48..57мм, статор60..100мм, одинаковые B–H/Br/токIамп5А/сеточные настройки/угловая сетка32положения. Изменена только дуга80→72градуса; объём PMуменьшаетсяна10%. Фактические сетки имеют разное число элементов вследствие изменения геометрии, локальные размеры одинаковы. Сталь линейнаяμ1000; потери и нагрев не рассчитываются. ЭДСосновнойгармоникиполучена проекцией потокосцепленийХХ на вращающийсябазис с циклическойсимметрией.','cases':results}
with open('magnet-arc-control-result.json' if N<32 else 'magnet-arc-result.json','w') as f:json.dump(out,f)
print('Magnet arc',[(r['full_magnet_arc_deg'],r['statistics']) for r in results])