Сферический магнитный экран: аналитическая проверка

"""Собственная учебная лабораторная FEMM 5: встроенный Python.
Команды FEMM доступны глобально, внешний пакет femm не используется.
Каждая сетка сохраняется в отдельный FEM, закрывается, открывается повторно
и решается с русскими именами свойств. Новый запуск создаёт calculated-NN.
Опубликованные модели и квитанции не изменяются. Размеры: мм; результаты: СИ.
"""
import os
import sys
import math
import json
CONFIG = {'a': 100, 'b': 110, 'B': 0.001, 'variants': [10, 100, 1000]}
MODE = 'sphere'
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 in ('wires','wireL','skin','proximity','busbars'):
    variants=CONFIG.get('variants',[0])
    for variant in variants:
      shapes=('solid','tube') if MODE=='busbars' else ('solid',)
      for shape in shapes:
       for level in range(3):
        f=variant if MODE in ('skin','proximity','busbars') else 0
        r=CONFIG.get('r',math.sqrt(20));d=CONFIG.get('d',30)
        if MODE=='skin':d=200
        outer=5000 if MODE=='wireL' else (1000 if MODE=='wires' else max(300,8*d))
        start(False,f,outer,level)
        if MODE in ('wires','wireL'):
            local=150 if MODE=='wires' else 750
            circle(0,0,local,4/(2**level))
            label(0,local*.7,'Воздух',(5 if MODE=='wires' else 15)/(2**level),3)
        sigma=CONFIG.get('sigma',0);material('Медь проводников',1,sigma)
        names=[]
        for side in (-1,1):
            name='Прямой проводник' if side<0 else 'Обратный проводник'
            mi_addcircprop(name,-side*CONFIG['I'],1);names.append(name)
            rr=6 if MODE=='busbars' and shape=='tube' else r
            circle(side*d/2,0,rr,4/(2**level))
            if MODE=='busbars' and shape=='tube':
                circle(side*d/2,0,4,4/(2**level));label(side*d/2,0,'Воздух',r/3/(2**level),2)
                label(side*d/2+5,0,'Медь проводников',r/4/(2**level),10,name,1)
            else:label(side*d/2,0,'Медь проводников',r/4/(2**level),10,name,1)
        label(0,outer*.6,'Воздух',outer/16/(2**level),1)
        row=solve(level,variant,'-'+shape);row['shape']=shape
        energy_identity(row,names,f)
        if MODE=='wires':
            x,y=-d/2,50;bx,by=magnetic_b(x,y);bx=real(bx);by=real(by)
            exactx=0;exacty=0
            for side in (-1,1):
                xx=x-side*d/2;ii=-side*CONFIG['I'];dist2=(xx*.001)**2+(y*.001)**2
                exactx-=MU0*ii*(y*.001)/(2*math.pi*dist2);exacty+=MU0*ii*(xx*.001)/(2*math.pi*dist2)
            row.update({'Bx_T':bx,'By_T':by,'analytic_Bx_T':exactx,'analytic_By_T':exacty})
            error=math.sqrt((bx-exactx)**2+(by-exacty)**2)/math.sqrt(exactx**2+exacty**2)
            check('Векторный закон Био — Савара '+row['model'],error,0.035/(2**level))
            row['primary']=math.sqrt(bx*bx+by*by)
        elif MODE=='wireL':
            inductance=2*row['energy_J']/CONFIG['I']**2
            exact=MU0/math.pi*(math.log(d/r)+.25)
            row.update({'L_H_per_m':inductance,'analytic_L_H_per_m':exact,'primary':inductance})
            check('Аналитическая индуктивность линии '+row['model'],abs(inductance/exact-1),0.03/(2**level))
        else:
            resistance=2*row['ohmic_loss_W']/CONFIG['I']**2
            row.update({'loop_R_ohm_per_m':resistance,'skin_only_loop_R_ohm_per_m':2*skin_resistance(r,sigma,f),'primary':resistance})
            if MODE=='skin':check('Внутреннее сопротивление по функциям Бесселя '+row['model'],abs(resistance/(2*skin_resistance(r,sigma,f))-1),0.018)
            elif MODE=='proximity':check('Потери не ниже двух изолированных проводов '+row['model'],max(0,2*skin_resistance(r,sigma,f)/resistance-1),0.015)
            else:
                dc=2/(math.pi*20e-6*sigma)
                check('AC-сопротивление не ниже DC '+row['model'],max(0,dc/resistance-1),0.01)
        records.append(row);finish()

elif MODE in ('helmholtz','maxwell','brooks','solenoid','mutual'):
    for variant in CONFIG['variants']:
      excitations=(1,2) if MODE=='mutual' else (1,)
      for excitation in excitations:
       for level in range(3):
        outer=2500 if MODE=='solenoid' else 1500
        start(True,0,outer,level);material('Медь обмотки')
        rect(0,-400,400,400)
        label(300,0,'Воздух',8/(2**level),3)
        names=[];specs=[]
        if MODE=='helmholtz':
            R=CONFIG['R'];sep=variant*R
            specs=[(R,-sep/2,2,2,CONFIG['I'],CONFIG['N']),(R,sep/2,2,2,CONFIG['I'],CONFIG['N'])]
        elif MODE=='maxwell':
            R=CONFIG['R'];rr=R*math.sqrt(4/7);zz=R*math.sqrt(3/7)*variant
            specs=[(R,0,2,2,CONFIG['I'],CONFIG['N']),(rr,-zz,2,2,CONFIG['I']*49/64,CONFIG['N']),(rr,zz,2,2,CONFIG['I']*49/64,CONFIG['N'])]
        elif MODE=='brooks':specs=[(1.5*variant,0,variant,variant,CONFIG['I'],CONFIG['N'])]
        elif MODE=='solenoid':specs=[(32.5,0,5,variant,CONFIG['I'],CONFIG['N'])]
        else:
            specs=[(50,-variant/2,4,4,1 if excitation==1 else 0,CONFIG['N']),(60,variant/2,4,4,1 if excitation==2 else 0,CONFIG['N'])]
        mesh=(variant/20 if MODE=='brooks' else 1)/(2**level)
        for num,spec in enumerate(specs,1):names.append(coil(*spec,mesh,num))
        label(outer*.5,0,'Воздух',outer/18/(2**level))
        row=solve(level,variant,'-e%d'%excitation);row['excitation']=excitation
        energy_identity(row,names,0)
        br,bz=magnetic_b(0,0);bz=real(bz);row['Bz_center_T']=bz
        if MODE in ('helmholtz','maxwell'):
            exact=sum(loop_field(r,z,n,i,0) for r,z,w,h,i,n in specs)
            row['analytic_Bz_center_T']=exact
            check('Знак и модуль осевого поля '+row['model'],abs(bz/exact-1),0.025/(2**level))
            profile=[]
            for frac in (-.3,-.2,-.1,0,.1,.2,.3):
                zz=frac*CONFIG['R'];bb=real(magnetic_b(0,zz)[1]);theory=sum(loop_field(r,z,n,i,zz) for r,z,w,h,i,n in specs)
                profile.append([zz,bb,theory]);check('Профиль поля z='+str(zz)+' '+row['model'],abs(bb/theory-1),0.035/(2**level))
            row['axis_profile_mm_T']=profile;row['uniformity_20percent']=max(abs(p[1]/bz-1) for p in profile if abs(p[0])<=.2*CONFIG['R']);row['primary']=bz
        elif MODE=='brooks':
            ll=2*row['energy_J']/CONFIG['I']**2;exact=1.6994e-6*1.5*variant*.001*CONFIG['N']**2
            row.update({'L_H':ll,'analytic_Grover_L_H':exact,'primary':ll})
            check('Формула Брукса — Гровера '+row['model'],abs(ll/exact-1),0.025)
        elif MODE=='solenoid':
            aa=variant*.001/2;exact=MU0*CONFIG['N']*CONFIG['I']/(2*.005)*(asinh(.035/aa)-asinh(.030/aa))
            ll=2*row['energy_J']/CONFIG['I']**2
            row.update({'analytic_center_T':exact,'L_H':ll,'primary':bz})
            check('Интеграл Био — Савара для конечной обмотки '+row['model'],abs(bz/exact-1),0.012)
        else:
            one=circuit_values(names[0]);two=circuit_values(names[1])
            ll=(one if excitation==1 else two)[2].real;mm=(two if excitation==1 else one)[2].real
            row.update({'L_self_H':ll,'M_H':mm,'primary':mm})
        records.append(row);finish()
    if MODE=='mutual':
        for variant in CONFIG['variants']:
            a=[r for r in records if r['variant']==variant and r['level']==3 and r['excitation']==1][0]
            b=[r for r in records if r['variant']==variant and r['level']==3 and r['excitation']==2][0]
            check('Взаимность M12=M21 при зазоре '+str(variant),abs(a['M_H']/b['M_H']-1),0.003)
            check('Положительная определённость матрицы L при зазоре '+str(variant),max(0,a['M_H']**2/(a['L_self_H']*b['L_self_H'])-1),1e-8)

elif MODE in ('cylinder','sphere','acshield'):
    for variant in CONFIG['variants']:
      for level in range(3):
        axis=MODE!='cylinder';a=CONFIG['a'];b=CONFIG['b'];mu=CONFIG.get('mu',variant);f=variant if MODE=='acshield' else 0
        sigma=CONFIG.get('sigma',0);outer=b*12
        start(axis,f,outer,level,CONFIG['B']);material('Материал экрана',mu,sigma)
        if axis:
            halfcircle(a,3/(2**level));halfcircle(b,3/(2**level))
            label((a+b)/2,0,'Материал экрана',(.5 if MODE=='acshield' else (b-a)/8)/(2**level),5)
        else:
            circle(0,0,a,3/(2**level));circle(0,0,b,3/(2**level));label((a+b)/2,0,'Материал экрана',(b-a)/6/(2**level),5)
        label(a/2,0,'Воздух',a/8/(2**level),2)
        label(b*3,0,'Воздух',outer/18/(2**level),1)
        row=solve(level,variant);bx,by=magnetic_b(0,0);bb=abs(complex(by if axis else bx))
        ratio=bb/CONFIG['B'];row.update({'B_inside_T':bb,'attenuation_ratio':ratio,'primary':ratio})
        if MODE=='cylinder':shield=((mu+1)**2-(mu-1)**2*(a/b)**2)/(4*mu)
        else:shield=((2*mu+1)*(mu+2)-2*(mu-1)**2*(a/b)**3)/(9*mu)
        row['analytic_DC_ratio']=1/shield
        if f==0:check('Аналитическое экранирование '+row['model'],abs(ratio*shield-1),0.018)
        if f>0:
            # Независимое интегрирование J^2/(2σ) по центрам элементов экрана.
            total=0
            for ei in range(1,mo_numelements()+1):
                elem=mo_getelement(ei)
                if int(elem[6])!=5:continue
                x,y,area=elem[3],elem[4],elem[5]
                values=mo_getpointvalues(x,y);jj=(complex(values[7])+complex(values[8]))*1e6
                volume=area*1e-6*(2*math.pi*x*.001 if axis else 1)
                total+=abs(jj)**2/(2*sigma)*volume
            row['element_Joule_integral_W']=total
            check('Интеграл Джоуля в экране '+row['model'],abs(total/row['ohmic_loss_W']-1),0.04)
        records.append(row);finish()
    if MODE=='acshield':
        last=[r for r in records if r['level']==3]
        for i in range(1,len(last)):
            check('Частотное усиление экранирования '+str(last[i]['variant']),max(0,last[i]['attenuation_ratio']/last[i-1]['attenuation_ratio']-1),0.01)

elif MODE in ('sensor','pipe','heating'):
    for variant in CONFIG['variants']:
      for level in range(3):
        f=variant if MODE=='heating' else CONFIG['f']
        start(True,f,500,level);material('Медь обмотки');material('Проводящая мишень',1,CONFIG['sigma'])
        if MODE=='sensor':
            name=coil(20,variant+3,4,6,CONFIG['I'],CONFIG['N'],1/(2**level),1)
            rect(0,-8,60,0);label(30,-4,'Проводящая мишень',1/(2**level),5)
        elif MODE=='heating':
            name=coil(40,10,6,8,CONFIG['I'],CONFIG['N'],1/(2**level),1)
            rect(0,-5,60,0);label(30,-2.5,'Проводящая мишень',.6/(2**level),5)
        else:
            name=coil(24,0,4,6,CONFIG['I'],CONFIG['N'],.8/(2**level),1)
            polygon([(16,-40),(20,-40),(20,-4),(20-variant,-4),(20-variant,4),(20,4),(20,40),(16,40)])
            label(18,-20,'Проводящая мишень',.8/(2**level),5)
        label(300,0,'Воздух',30/(2**level))
        row=solve(level,variant);energy_identity(row,[name],f)
        ii,vv,pp=circuit_values(name);zz=vv/ii
        row.update({'impedance_real_ohm':zz.real,'impedance_imag_ohm':zz.imag,'L_apparent_H':zz.imag/(2*math.pi*f),'primary':row['ohmic_loss_W']})
        check('Пассивность мишени '+row['model'],max(0,-zz.real),1e-10)
        check('Реактивное сопротивление катушки '+row['model'],max(0,-zz.imag),1e-10)
        records.append(row);finish()
    finest=[r for r in records if r['level']==3]
    if MODE=='sensor':
        for i in range(1,len(finest)):check('Убывание потерь с зазором '+str(finest[i]['variant']),max(0,finest[i]['ohmic_loss_W']/finest[i-1]['ohmic_loss_W']-1),0.005)
    if MODE=='pipe':
        difference=abs(finest[1]['impedance_imag_ohm']/finest[0]['impedance_imag_ohm']-1)
        check('Утонение даёт разрешимый сигнал',max(0,0.001-difference),0)
        records.append({'measurement':'thinning_signal','relative_impedance_change':difference})

# Сходимость по трём сеткам проверяется отдельно для каждого физического режима.
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 ('acshield','heating','pipe','sensor','proximity','skin','busbars') else 0.012)
result={'runtime':'FEMM 5; embedded MicroPython','case_id':'magnetic-shield-sphere','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']