Комплексная матрица импедансов: токи, общая оболочка и нагрузки

"""Собственная учебная лабораторная FEMM 5: встроенный Python.
Команды FEMM доступны глобально, внешний пакет femm не используется.
Каждая сетка сохраняется в отдельный FEM, закрывается, открывается повторно
и решается с русскими именами свойств. Новый запуск создаёт calculated-NN.
Опубликованные модели и квитанции не изменяются. Размеры: мм; результаты: СИ.
"""
import os
import sys
import math
import json
CONFIG = {'sigma': 58000000.0, 'variants': [50, 5000, 50000], 'loads': [[0.005, 0.05], [0.05, 0.005]]}
MODE = 'matrix'
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)))


def pair(z):return [z.real,z.imag]
def phase(z):return math.atan2(z.imag,z.real)*180/math.pi
def group_current(group):
    mo_groupselectblock(group);current=complex(mo_blockintegral(7));mo_clearblock();return current
def circuit_power(names):
    total=0
    for name in names:
        ii,vv,pp=circuit_values(name);total+=(vv*ii.conjugate()).real/2
    return total

if MODE=='shade':
 def build_shade(f,level,currents):
    newdocument(0);mi_probdef(f,'millimeters','planar',1e-10,100,30,0)
    material('Воздух');material('Линейный магнитопровод',500);material('Изолированная обмотка');material('Медь экранирующего витка',1,CONFIG['sigma'])
    mi_addboundprop('Удалённая граница A=0',0,0,0,0,0,0,0,0,0)
    rect(-40,-15,40,60)
    for x,y in ((0,-15),(40,20),(0,60),(-40,20)):
        mi_selectsegment(x,y);mi_setsegmentprop('Удалённая граница A=0',4/(2**level),0,0,0);mi_clearselected()
    polygon([(-18,0),(18,0),(18,36),(12,36),(12,6),(-12,6),(-12,36),(-18,36)])
    label(0,3,'Линейный магнитопровод',.8/(2**level),5)
    rect(-18,38,18,44);label(0,41,'Линейный магнитопровод',.8/(2**level),6)
    names=['Первичная обмотка','Экранирующий виток']
    for name,current in zip(names,currents):mi_addcircprop(name,current,1)
    for x,sign in ((-20,1),(-10,-1)):
        rect(x-1,12,x+1,24);label(x,18,'Изолированная обмотка',.5/(2**level),11,names[0],200*sign)
    for x,sign in ((-19.5,1),(-15,-1)):
        rect(x-.5,30,x+.5,35);label(x,32.5,'Медь экранирующего витка',.25/(2**level),12,names[1],sign)
    # Прорезь выходит на поверхность полюса: без неё железная перемычка
    # вокруг внутренней стороны витка шунтировала бы его полезную связь.
    rect(-15.5,35,-14.5,36);label(-15,35.5,'Воздух',.2/(2**level),3)
    # Дополнительные линии задают подробную сетку зазора и полюсных краёв.
    rect(-24,27,24,48);label(22,37,'Воздух',.5/(2**level),3)
    label(-17,29,'Линейный магнитопровод',.4/(2**level),5);label(15,30,'Линейный магнитопровод',.4/(2**level),5)
    label(0,20,'Воздух',1.5/(2**level),2);label(30,0,'Воздух',4/(2**level),1)
    return names
 def pole_average(x1,x2):
    aa=complex(mo_getpointvalues(x1,37)[0]);bb=complex(mo_getpointvalues(x2,37)[0])
    return (aa-bb)/((x2-x1)*.001)
 for variant in CONFIG['variants']:
  for level in range(3):
    Z=[[0j,0j],[0j,0j]]
    for excitation in (0,1):
        currents=[1 if j==excitation else 0 for j in range(2)]
        names=build_shade(variant,level,currents)
        row=solve(level,variant,'-exc%d'%excitation);row['excitation']=excitation;row['shape']='matrix'
        energy_identity(row,names,variant)
        for j,name in enumerate(names):Z[j][excitation]=circuit_values(name)[1]
        row.update({'matrix_column_complex_ohm':[pair(Z[j][excitation]) for j in range(2)],'primary':abs(Z[excitation][excitation])})
        if excitation==0:
            row['open_shaded_B_complex_T']=pair(pole_average(-18,-15.5));row['open_unshaded_B_complex_T']=pair(pole_average(-14.5,-12))
        records.append(row);finish()
    check('Взаимность матрицы экранирующего витка '+str(variant)+' сетка '+str(level+1),abs(Z[0][1]-Z[1][0])/max(abs(Z[0][1]),abs(Z[1][0])),.006)
    current=-Z[1][0]/(Z[1][1]+CONFIG['Rend'])
    names=build_shade(variant,level,[1,current])
    row=solve(level,variant,'-shorted');row['excitation']=2;row['shape']='shorted'
    energy_identity(row,names,variant)
    ii,vv,pp=circuit_values(names[1]);primary=circuit_values(names[0]);bend=CONFIG['Rend']*abs(current)**2/2
    shaded=pole_average(-18,-15.5);unshaded=pole_average(-14.5,-12)
    delta=phase(shaded/unshaded)
    kvl=abs(vv+CONFIG['Rend']*ii)/max(abs(Z[1][0]),1e-20)
    power=(primary[1]*primary[0].conjugate()).real/2
    check('Кирхгоф замкнутого витка '+row['model'],kvl,.0002)
    check('Мощность источника равна потерям поля и торцов '+row['model'],abs(power-row['ohmic_loss_W']-bend)/max(abs(power),1e-20),.015)
    check('Пассивность короткозамкнутого витка '+row['model'],max(0,-power),0)
    row.update({'short_turn_current_real_A':current.real,'short_turn_current_imag_A':current.imag,'short_turn_current_amplitude_A':abs(current),'shaded_field_amplitude_T':abs(shaded),'unshaded_field_amplitude_T':abs(unshaded),'spatial_phase_shift_deg':delta,'end_return_loss_W':bend,'source_active_power_W':power,'loop_KVL_relative':kvl,'primary':abs(current)})
    records.append(row);finish()
 finest=[r for r in records if r.get('shape')=='shorted' and r['level']==3]
 for row in finest:
    check('Разрешимый пространственный фазовый сдвиг '+str(row['variant']),max(0,1-abs(row['spatial_phase_shift_deg'])),0)

elif MODE=='matrix':
 def build_matrix(f,level,currents):
    start(False,f,40,level);material('Медь жил и общей оболочки',1,CONFIG['sigma'])
    names=[]
    allcurrents=currents+[-sum(currents)]
    for j,current in enumerate(allcurrents):
        name='Жила %d'%(j+1) if j<3 else 'Общая возвратная оболочка'
        mi_addcircprop(name,current,1);names.append(name)
    for j,x in enumerate((-4,0,4)):
        circle(x,0,1.5,.8/(2**level));label(x,0,'Медь жил и общей оболочки',.24/(2**level),11+j,names[j])
    circle(0,0,10,1.2/(2**level));circle(0,0,11,1.2/(2**level));label(0,10.5,'Медь жил и общей оболочки',.2/(2**level),14,names[3])
    label(0,5,'Воздух',.7/(2**level),3);label(0,25,'Воздух',4/(2**level),1)
    return names
 def get_ports(names):
    values=[circuit_values(name) for name in names]
    return [values[j][1]-values[3][1] for j in range(3)]
 for variant in CONFIG['variants']:
  for level in range(3):
    Z=[[0j]*3 for j in range(3)]
    for excitation in range(3):
        currents=[1 if j==excitation else 0 for j in range(3)]
        names=build_matrix(variant,level,currents)
        row=solve(level,variant,'-exc%d'%excitation);row['excitation']=excitation;row['shape']='matrix'
        energy_identity(row,names,variant)
        ports=get_ports(names)
        for j in range(3):Z[j][excitation]=ports[j]
        for j,wanted in enumerate(currents+[-sum(currents)]):
            check('Интеграл плотности тока проводника '+str(j+1)+' '+row['model'],abs(group_current(11+j)-wanted),.002)
        dc=1/(CONFIG['sigma']*math.pi*(.0015)**2)+1/(CONFIG['sigma']*math.pi*(.011**2-.010**2))
        check('Нижняя граница равномерного DC-тока '+row['model'],max(0,dc/Z[excitation][excitation].real-1),.002)
        row.update({'port_impedances_complex_ohm':[pair(v) for v in ports],'port_self_resistance_ohm':Z[excitation][excitation].real,'port_self_reactance_ohm':Z[excitation][excitation].imag,'dc_uniform_lower_bound_ohm':dc,'primary':abs(ports[excitation])})
        records.append(row);finish()
    scale=max(abs(z) for line in Z for z in line)
    for j in range(3):
        for k in range(j):check('Комплексная взаимность Z%d%d '%(j+1,k+1)+str(variant)+' сетка '+str(level+1),abs(Z[j][k]-Z[k][j])/scale,.002)
    # Пассивность проверяется для независимого комплексного набора токов.
    probe=[1,complex(-.4,.8),complex(.2,-.5)]
    p=sum((sum(Z[j][k]*probe[k] for k in range(3))*probe[j].conjugate()).real/2 for j in range(3))
    check('Пассивность квадратичной формы '+str(variant)+' сетка '+str(level+1),max(0,-p),0)
    for load_index,loads in enumerate(CONFIG['loads']):
        aa=Z[1][1]+loads[0];bb=Z[1][2];cc=Z[2][1];dd=Z[2][2]+loads[1]
        det=aa*dd-bb*cc
        currents=[1,(-Z[1][0]*dd+bb*Z[2][0])/det,(-aa*Z[2][0]+cc*Z[1][0])/det]
        names=build_matrix(variant,level,currents)
        row=solve(level,variant,'-load%d'%load_index);row['excitation']=3+load_index;row['shape']='load%d'%load_index
        energy_identity(row,names,variant)
        ports=get_ports(names);predicted=[sum(Z[j][k]*currents[k] for k in range(3)) for j in range(3)]
        mismatch=max(abs(ports[j]-predicted[j]) for j in range(3))/max(abs(z) for z in ports)
        kvl=max(abs(ports[j]+loads[j-1]*currents[j]) for j in (1,2))/max(abs(z) for z in ports)
        loadloss=sum(loads[j-1]*abs(currents[j])**2/2 for j in (1,2));power=ports[0].real/2
        check('Суперпозиция полной комплексной матрицы '+row['model'],mismatch,.0002)
        check('Кирхгоф двух нагруженных портов '+row['model'],kvl,.0002)
        check('Мощность источника — металл и нагрузки '+row['model'],abs(power-row['ohmic_loss_W']-loadloss)/max(abs(power),1e-20),.015)
        row.update({'load_index':load_index,'load_resistances_ohm':loads,'port_currents_complex_A':[pair(z) for z in currents],'port_current_amplitudes_A':[abs(z) for z in currents],'port_voltages_complex_V':[pair(z) for z in ports],'load_losses_W':loadloss,'source_active_power_W':power,'matrix_superposition_relative':mismatch,'loaded_KVL_relative':kvl,'primary':power})
        records.append(row);finish()
    if level==2:records.append({'measurement':'full_complex_port_impedance_matrix','frequency_Hz':variant,'matrix_complex_ohm':[[pair(z) for z in line] for line in Z]})
# Сходимость по трём сеткам проверяется отдельно для каждого физического режима.
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 ('shade','matrix') else 0.012)
result={'runtime':'FEMM 5; embedded MicroPython','case_id':'ac-impedance-matrix-loads','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']