Сила на малой ферромагнитной частице в поле катушки

"""Собственная учебная лабораторная FEMM 5: встроенный Python.
Команды FEMM доступны глобально, внешний пакет femm не используется.
Каждая сетка сохраняется в отдельный FEM, закрывается, открывается повторно
и решается с русскими именами свойств. Новый запуск создаёт calculated-NN.
Опубликованные модели и квитанции не изменяются. Размеры: мм; результаты: СИ.
"""
import os
import sys
import math
import json
CONFIG = {'I': 1, 'N': 100, 'mu': 100, 'variants': [10, 20, 30]}
MODE = 'particle'
saved=globals().get('__file__','') or (sys.argv[0] if sys.argv else '')
base=os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd()
index=1
while os.path.exists(base+'/calculated-%02d'%index):index+=1
directory=base+'/calculated-%02d'%index
os.mkdir(directory)
records=[]
checks=[]
MU0=4*math.pi*1e-7

def real(x):
    if isinstance(x,(list,tuple)):x=x[0]
    return complex(x).real

def check(name,value,limit):
    checks.append({'name':name,'value':value,'limit':limit,'passed':math.isfinite(value) and value<=limit})

def material(name,mu=1,sigma=0):
    mi_addmaterial(name,mu,mu,0,0,sigma/1e6)

def label(x,y,name,mesh,group=1,circuit='<None>',turns=1):
    mi_addblocklabel(x,y);mi_selectlabel(x,y)
    mi_setblockprop(name,0,mesh,circuit,0,group,turns);mi_clearselected()

def polygon(points,boundary=None,spacing=0):
    for x,y in points:mi_addnode(x,y)
    for i,(x,y) in enumerate(points):
        x2,y2=points[(i+1)%len(points)];mi_addsegment(x,y,x2,y2)
        if boundary:
            mi_selectsegment((x+x2)/2,(y+y2)/2)
            mi_setsegmentprop(boundary,spacing,0,0,0);mi_clearselected()

def rect(x1,y1,x2,y2):polygon([(x1,y1),(x2,y1),(x2,y2),(x1,y2)])

def circle(x,y,r,angle,boundary=None):
    mi_addnode(x-r,y);mi_addnode(x+r,y)
    mi_addarc(x-r,y,x+r,y,180,angle);mi_addarc(x+r,y,x-r,y,180,angle)
    if boundary:
        mi_selectarcsegment(x,y-r);mi_selectarcsegment(x,y+r)
        mi_setarcsegmentprop(angle,boundary,0,0);mi_clearselected()

def halfcircle(r,angle,boundary=None):
    mi_addnode(0,-r);mi_addnode(0,r)
    mi_addarc(0,-r,0,r,180,angle)
    if boundary:
        mi_selectarcsegment(r,0);mi_setarcsegmentprop(angle,boundary,0,0);mi_clearselected()

def start(axis,f,outer,level,uniform=0):
    newdocument(0)
    mi_probdef(f,'millimeters','axi' if axis else 'planar',1e-10,1000,30,0)
    material('Воздух')
    # Коэффициенты пространственной функции A вводятся на миллиметр геометрии.
    mi_addboundprop('Внешняя граница',0,uniform/2000 if axis else 0,0 if axis else uniform/1000,0,0,0,0,0,0)
    if axis:
        # Ось r=0, z снизу вверх: исходные координаты библиотек не используются.
        halfcircle(outer,4/(2**level),'Внешняя граница')
        mi_addsegment(0,-outer,0,outer)
    else:circle(0,0,outer,4/(2**level),'Внешняя граница')

def close_axis():
    # FEMM автоматически разбивает осевой отрезок вершинами внутренних дуг.
    pass

def coil(r,z,width,height,current,turns,mesh,num):
    name='Обмотка %d'%num
    mi_addcircprop(name,current,1)
    rect(r-width/2,z-height/2,r+width/2,z+height/2)
    label(r,z,'Медь обмотки',mesh,10+num,name,turns)
    return name

def solve(level,variant,tag=''):
    name='model-%s-%02d%s.fem'%(str(variant).replace('.','p'),level+1,tag)
    mi_saveas(directory+'/'+name);mi_close();opendocument(directory+'/'+name)
    mi_analyze(1);mi_loadsolution();mo_smooth('off')
    mo_groupselectblock();energy=real(mo_blockintegral(2));loss=real(mo_blockintegral(4));mo_clearblock()
    return {'model':name,'level':level+1,'variant':variant,'nodes':mo_numnodes(),'elements':mo_numelements(),'energy_J':energy,'ohmic_loss_W':loss,'saved_reopened_solved':True}

def finish():mo_close();mi_close()

def loop_field(r,z,turns,current,position):
    r*=0.001;dz=(position-z)*0.001
    return MU0*turns*current*r*r/(2*(r*r+dz*dz)**1.5)

def bessel(n,z):
    term=(z/2)**n
    total=term
    for k in range(1,100):
        term*=-(z*z/4)/(k*(k+n));total+=term
        if abs(term)<1e-15*max(1,abs(total)):break
    return total

def skin_resistance(radius,sigma,f):
    a=radius*.001
    if f==0:return 1/(math.pi*a*a*sigma)
    delta=math.sqrt(2/(2*math.pi*f*MU0*sigma));k=(1-1j)/delta
    return (k*bessel(0,k*a)/(2*math.pi*a*sigma*bessel(1,k*a))).real

def circuit_values(name):
    i,v,psi=mo_getcircuitproperties(name)
    return complex(i),complex(v),complex(psi)

def magnetic_b(x,y):
    mo_smooth('on')
    values=mo_getpointvalues(x,y)
    mo_smooth('off')
    return values[1],values[2]

def asinh(x):return math.log(x+math.sqrt(x*x+1))

def energy_identity(row,circuits,frequency):
    power=0;fluxenergy=0
    for name in circuits:
        i,v,psi=circuit_values(name)
        power+=(v*i.conjugate()).real/2
        fluxenergy+=(psi*i.conjugate()).real/(4 if frequency else 2)
    row['circuit_active_power_W']=power
    row['circuit_energy_J']=fluxenergy
    if frequency:
        error=abs(power-row['ohmic_loss_W'])/max(abs(power),1e-20)
        row['power_balance_relative']=error
        check('Баланс активной мощности '+row['model'],error,0.015)
    else:
        error=abs(fluxenergy-row['energy_J'])/max(abs(fluxenergy),1e-20)
        row['energy_identity_relative']=error
        check('Энергия и потокосцепление '+row['model'],error,0.012/(2**(row['level']-1)))


if MODE=='slot':
 for variant in CONFIG['variants']:
  for enclosed in (0,1):
   for level in range(3):
    start(False,variant,300,level);material('Медь проводника',1,CONFIG['sigma']);material('Материал паза',1000 if enclosed else 1,0)
    polygon([(-20,-25),(20,-25),(20,25),(8,25),(8,-10),(-8,-10),(-8,25),(-20,25)])
    label(0,-18,'Материал паза',2/(2**level),5)
    names=[]
    for num,(x,current) in enumerate(((0,CONFIG['I']),(70,-CONFIG['I'])),1):
        name='Цепь проводника %d'%num;mi_addcircprop(name,current,1);names.append(name)
        rect(x-3,-6,x+3,6);label(x,0,'Медь проводника',.8/(2**level),10+num,name)
    circle(0,0,110,3/(2**level));label(0,70,'Воздух',3/(2**level),3);label(0,220,'Воздух',25/(2**level))
    row=solve(level,variant,'-core%d'%enclosed);row['core']=enclosed
    energy_identity(row,names,variant);row['R_loop_ohm']=2*row['ohmic_loss_W']/CONFIG['I']**2;row['primary']=row['R_loop_ohm']
    dc=2/(CONFIG['sigma']*.006*.012)
    check('Нижняя DC-граница '+row['model'],max(0,dc/row['R_loop_ohm']-1),.01)
    row['shape']=str(enclosed);records.append(row);finish()

elif MODE=='lamination':
 for variant in CONFIG['variants']:
  for level in range(3):
    f=CONFIG['f'];thickness=variant;width=20;mesh=thickness/5/(2**level)
    newdocument(0);mi_probdef(f,'millimeters','planar',1e-11,1000,30,0)
    material('Электротехническая пластина',CONFIG['mu'],CONFIG['sigma'])
    aa=CONFIG['B']*thickness*.001/2
    mi_addboundprop('Потенциал верхней поверхности',aa,0,0,0,0,0,0,0,0)
    mi_addboundprop('Потенциал нижней поверхности',-aa,0,0,0,0,0,0,0,0)
    # Осевой повторяющийся фрагмент с нулевым потоком через боковые торцы.
    rect(-width/2,-thickness/2,width/2,thickness/2)
    for y,name in ((thickness/2,'Потенциал верхней поверхности'),(-thickness/2,'Потенциал нижней поверхности')):
        mi_selectsegment(0,y);mi_setsegmentprop(name,1,0,0,0);mi_clearselected()
    label(0,0,'Электротехническая пластина',mesh,5)
    row=solve(level,variant)
    volume=width*thickness*1e-6
    thin=math.pi**2/6*CONFIG['sigma']*f*f*CONFIG['B']**2*(thickness*.001)**2*volume
    row.update({'thin_sheet_loss_W':thin,'volume_m3':volume,'loss_density_W_m3':row['ohmic_loss_W']/volume,'primary':row['ohmic_loss_W']})
    check('Тонколистовое аналитическое решение '+row['model'],abs(row['ohmic_loss_W']/thin-1),.025)
    records.append(row);finish()

elif MODE=='litz':
 for variant in CONFIG['variants']:
  for balanced in (0,1):
   for level in range(3):
    start(False,variant,50,level);material('Медь изолированных жил',1,CONFIG['sigma'])
    positions=[(0,0)]+[(math.cos(k*math.pi/3),math.sin(k*math.pi/3)) for k in range(6)]
    names=[]
    if not balanced:
        mi_addcircprop('Общие выводы пучка',CONFIG['I'],0);names.append('Общие выводы пучка')
    for num,(x,y) in enumerate(positions):
        name='Ток жилы %d'%(num+1) if balanced else 'Общие выводы пучка'
        if balanced:mi_addcircprop(name,CONFIG['I']/7,1);names.append(name)
        circle(x,y,.4,4/(2**level));label(x,y,'Медь изолированных жил',.1/(2**level),10+num,name)
    mi_addcircprop('Обратный провод',-CONFIG['I'],1);names.append('Обратный провод')
    circle(10,0,.4*math.sqrt(7),4/(2**level));label(10,0,'Медь изолированных жил',.15/(2**level),30,'Обратный провод')
    circle(0,0,15,4/(2**level));label(0,7,'Воздух',.5/(2**level),3);label(0,35,'Воздух',5/(2**level))
    row=solve(level,variant,'-balanced%d'%balanced);row['balanced']=balanced
    energy_identity(row,names,variant)
    currents=[]
    for num in range(7):
        mo_groupselectblock(10+num);currents.append(complex(mo_blockintegral(7)));mo_clearblock()
    check('Суммарный ток семи жил '+row['model'],abs(sum(currents)-CONFIG['I'])/CONFIG['I'],1e-6)
    row['strand_current_amplitudes_A']=[abs(i) for i in currents]
    row['current_spread_relative']=(max(abs(i) for i in currents)-min(abs(i) for i in currents))/(CONFIG['I']/7)
    row['primary']=row['ohmic_loss_W'];row['shape']=str(balanced)
    records.append(row);finish()

elif MODE=='threephase':
 for variant in CONFIG['variants']:
  for shield in (0,1):
   for level in range(3):
    start(False,variant,350,level);material('Медь фазных шин',1,CONFIG['sigma']);material('Материал экрана',1,35e6 if shield else 0)
    circle(0,0,75,3/(2**level));circle(0,0,80,3/(2**level));label(77,0,'Материал экрана',1/(2**level),5)
    names=[];current_sum=0j
    for num,(x,phase) in enumerate(((-20,0),(0,-120),(20,120))):
        name=('A','B','C')[num]+' — фазная шина';current=CONFIG['I']*complex(math.cos(phase*math.pi/180),math.sin(phase*math.pi/180))
        mi_addcircprop(name,current,1);names.append(name);current_sum+=current
        rect(x-3,-7.5,x+3,7.5);label(x,0,'Медь фазных шин',1/(2**level),10+num,name)
    if shield:
        mi_addcircprop('Экран без внешнего тока',0,1);mi_selectlabel(77,0);mi_setblockprop('Материал экрана',0,1/(2**level),'Экран без внешнего тока',0,5,1);mi_clearselected();names.append('Экран без внешнего тока')
    label(0,40,'Воздух',3/(2**level),3);label(0,200,'Воздух',30/(2**level))
    row=solve(level,variant,'-shield%d'%shield);row['shield']=shield;row['shape']=str(shield)
    energy_identity(row,names,variant);check('Баланс фазных токов '+row['model'],abs(current_sum)/CONFIG['I'],1e-12)
    mo_groupselectblock(5);row['shield_loss_W']=real(mo_blockintegral(4));mo_clearblock()
    row['primary']=row['ohmic_loss_W'];records.append(row);finish()

elif MODE=='line-matrix':
 for variant in CONFIG['variants']:
  for level in range(3):
    start(False,0,800,level);material('Немагнитные провода')
    names=[]
    for num,x in enumerate((-30,-10,10,30)):
        name=('A','B','C','N')[num]+' — провод линии'
        current=1 if num==variant else (-1 if num==3 else 0)
        mi_addcircprop(name,current,1);names.append(name)
        circle(x,0,3,3/(2**level));label(x,0,'Немагнитные провода',.5/(2**level),10+num,name)
    circle(0,0,100,3/(2**level));label(0,70,'Воздух',3/(2**level),3);label(0,600,'Воздух',50/(2**level))
    row=solve(level,variant);energy_identity(row,names,0)
    flux=[circuit_values(n)[2].real for n in names]
    column=[flux[k]-flux[3] for k in range(3)];row['matrix_column_H_per_m']=column
    row['self_L_H_per_m']=column[variant];row['primary']=column[variant]
    distance=abs((-30,-10,10)[variant]-30);exact=MU0/math.pi*(math.log(distance/3)+.25)
    row['analytic_self_H_per_m']=exact
    check('Аналитическая петля с нейтралью '+row['model'],abs(column[variant]/exact-1),.025)
    records.append(row);finish()
 matrix=[[0,0,0] for _ in range(3)]
 for row in records:
    if row['level']==3:
        for k in range(3):matrix[k][row['variant']]=row['matrix_column_H_per_m'][k]
 for i in range(3):
  for j in range(3):check('Взаимность матрицы L%d%d'%(i,j),abs(matrix[i][j]-matrix[j][i])/math.sqrt(matrix[i][i]*matrix[j][j]),.002)
 for i in range(3):
  for j in range(i):check('Положительность главного минора L%d%d'%(i,j),max(0,matrix[i][j]**2/(matrix[i][i]*matrix[j][j])-1),1e-8)
 average=[[sum(matrix[(i+k)%3][(j+k)%3] for k in range(3))/3 for j in range(3)] for i in range(3)]
 records.append({'measurement':'matrix_and_ideal_transposition','matrix_H_per_m':matrix,'transposed_average_H_per_m':average})

elif MODE=='earth':
 delta=math.sqrt(2/(2*math.pi*CONFIG['f']*MU0*CONFIG['sigma']))
 for variant in CONFIG['variants']:
  for level in range(3):
    outer=delta*variant*1000
    start(False,CONFIG['f'],outer,level);material('Медь надземного провода',1,57e6);material('Проводящая среда возврата',1,CONFIG['sigma'])
    mi_addnode(-outer,0);mi_addnode(outer,0);mi_addsegment(-outer,0,outer,0)
    # Две смежные области грунта — один проводник с ОБЩИМ полным током.
    # Параллельный тип делит ток по областям вместо повторения -I в каждой.
    mi_addcircprop('Надземный провод',CONFIG['I'],1);mi_addcircprop('Возврат через среду',-CONFIG['I'],0)
    circle(0,2000,20,4/(2**level));label(0,2000,'Медь надземного провода',2/(2**level),10,'Надземный провод')
    circle(0,0,20000,4/(2**level))
    label(0,10000,'Воздух',500/(2**level),3)
    label(0,-10000,'Проводящая среда возврата',1000/(2**level),5,'Возврат через среду')
    label(0,outer*.7,'Воздух',outer/16/(2**level),1)
    label(0,-outer*.7,'Проводящая среда возврата',delta*1000/8/(2**level),6,'Возврат через среду')
    row=solve(level,variant);energy_identity(row,['Надземный провод','Возврат через среду'],CONFIG['f'])
    i,v,p=circuit_values('Надземный провод');ig,vg,pg=circuit_values('Возврат через среду')
    zz=(v-vg)/i
    row.update({'skin_depth_m':delta,'outer_radius_m':outer*.001,'loop_Zreal_ohm_per_m':zz.real,'loop_Zimag_ohm_per_m':zz.imag,'primary':zz.real})
    ground_current=0j
    for group in (5,6):
        mo_groupselectblock(group);ground_current+=complex(mo_blockintegral(7));mo_clearblock()
    row['integrated_ground_current_real_A']=ground_current.real
    row['integrated_ground_current_imag_A']=ground_current.imag
    check('Интегральный баланс токов провод — среда '+row['model'],abs(i+ground_current)/CONFIG['I'],1e-6)
    records.append(row);finish()
 finest=[r for r in records if r['level']==3]
 check('Удаление внешней границы от6δдо8δ',abs(finest[-1]['primary']/finest[-2]['primary']-1),.03)

elif MODE=='pipe-joint':
 for variant in CONFIG['variants']:
  for level in range(3):
    start(True,CONFIG['f'],600,level);material('Медь обмотки');material('Основной металл трубы',1,6e6);material('Кольцевая зона трубы',1,6e6*variant)
    name=coil(48,0,4,20,CONFIG['I'],CONFIG['N'],1/(2**level),1)
    for lo,hi,mat,group in ((-40,-4,'Основной металл трубы',5),(-4,4,'Кольцевая зона трубы',6),(4,40,'Основной металл трубы',5)):
        rect(35,lo,40,hi);label(37,(lo+hi)/2,mat,.5/(2**level),group)
    label(0.5,0,'Воздух',2/(2**level),3);label(300,0,'Воздух',30/(2**level),1)
    row=solve(level,variant);energy_identity(row,[name],CONFIG['f'])
    for group,key in ((5,'bulk_loss_W'),(6,'joint_loss_W')):
        mo_groupselectblock(group);row[key]=real(mo_blockintegral(4));mo_clearblock()
    check('Сумма потерь зон трубы '+row['model'],abs((row['bulk_loss_W']+row['joint_loss_W'])/row['ohmic_loss_W']-1),1e-8)
    row['joint_loss_density_W_m3']=row['joint_loss_W']/(math.pi*(.04**2-.035**2)*.008)
    row['primary']=row['joint_loss_W'];records.append(row);finish()

elif MODE=='particle':
 for variant in CONFIG['variants']:
  for level in range(3):
    start(True,0,400,level);material('Медь обмотки');material('Ферромагнитная частица',CONFIG['mu'])
    name=coil(20,0,2,2,CONFIG['I'],CONFIG['N'],.5/(2**level),1)
    for radius in (1,3):
        mi_addnode(0,variant-radius);mi_addnode(0,variant+radius);mi_addarc(0,variant-radius,0,variant+radius,180,2/(2**level))
    label(.5,variant,'Ферромагнитная частица',.2/(2**level),5)
    label(2,variant,'Воздух',.2/(2**level),6)
    rect(0,-50,60,80);label(40,0,'Воздух',1.5/(2**level),3);label(200,0,'Воздух',25/(2**level))
    row=solve(level,variant);energy_identity(row,[name],0)
    mo_groupselectblock(5);force=real(mo_blockintegral(19));volume=real(mo_blockintegral(10));mo_clearblock()
    r=.02;z=variant*.001;b=MU0*CONFIG['N']*CONFIG['I']*r*r/(2*(r*r+z*z)**1.5);deriv=-3*z/(r*r+z*z)*b
    exact=3*(4/3*math.pi*.001**3)*(CONFIG['mu']-1)/(CONFIG['mu']+2)/MU0*b*deriv
    row.update({'Fz_N':force,'dipole_force_N':exact,'particle_volume_m3':volume,'primary':force})
    check('Дипольная оценка малой сферы '+row['model'],abs(force/exact-1),.04)
    check('Сила направлена к катушке '+row['model'],max(0,force),0)
    records.append(row);finish()
# Сходимость по трём сеткам проверяется отдельно для каждого физического режима.
groups={}
for row in records:
    if 'level' not in row:continue
    key=str(row['variant'])+str(row.get('shape',''))+str(row.get('excitation',''))
    if key not in groups:groups[key]=[]
    groups[key].append(row)
for key,rows in groups.items():
    rows.sort(key=lambda row:row['level'])
    value=abs(rows[-1]['primary']/rows[-2]['primary']-1)
    check('Сходимость измеряемой величины, режим '+key,value,0.025 if MODE in ('slot','lamination','litz','threephase','line-matrix','earth','pipe-joint','particle') else 0.012)
result={'runtime':'FEMM 5; embedded MicroPython','case_id':'ferro-particle-force','records':records,'checks':checks,'passed':all(c['passed'] for c in checks),'scope':'Самостоятельная учебная 2D/осесимметричная модель; численные проверки относятся к указанным параметрам.'}
with open(directory+'/result.json','w') as stream:json.dump(result,stream)
print(json.dumps(result))
assert result['passed']