Резонансная передача энергии: реальные витки, компенсация и КПД

"""Собственная учебная лабораторная FEMM 5: встроенный Python.
Команды FEMM доступны глобально, внешний пакет femm не используется.
Каждая сетка сохраняется в отдельный FEM, закрывается, открывается повторно
и решается с русскими именами свойств. Новый запуск создаёт calculated-NN.
Опубликованные модели и квитанции не изменяются. Размеры: мм; результаты: СИ.
"""
import os
import sys
import math
import json
CONFIG = {}
MODE = 'resonance'
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 build(f,level,currents):
    start(True,f,250,level);material('Медь отдельных изолированных витков',1,58e6)
    names=['Передающая катушка','Приёмная катушка']
    for name,current in zip(names,currents):mi_addcircprop(name,current,1)
    for winding,z in enumerate((-10,10)):
        for turn in range(10):
            radius=20+2*turn
            circle(radius,z,.5,3/(2**level));label(radius,z,'Медь отдельных изолированных витков',.15/(2**level),11+winding,names[winding],1)
    rect(0,-20,50,20);label(45,0,'Воздух',1.5/(2**level),3);label(180,0,'Воздух',20/(2**level),1)
    return names

RLEAD=.01
RLOAD=2
VOLTAGE=1
F0=100000
calibration=[]
for level in range(3):
    names=build(F0,level,[1,0]);row=solve(level,F0,'-calibration');row['shape']='calibration';row['excitation']=0
    energy_identity(row,names,F0)
    z=circuit_values(names[0])[1];inductance=z.imag/(2*math.pi*F0)
    row.update({'inductance_calibration_H':inductance,'primary':inductance});records.append(row);calibration.append(inductance);finish()
C=1/((2*math.pi*F0)**2*calibration[-1])
dc=sum(2*math.pi*(.020+.002*k)/(58e6*math.pi*.0005**2) for k in range(10))
cache={}
def measure(f,level,refinement=False):
    if (f,level) in cache:return cache[(f,level)]
    Z=[[0j,0j],[0j,0j]]
    for excitation in (0,1):
        names=build(f,level,[1 if k==excitation else 0 for k in range(2)])
        row=solve(level,f,'-exc%d'%excitation);row['shape']='refinement-matrix' if refinement else 'matrix';row['excitation']=excitation
        energy_identity(row,names,f)
        for k in range(2):Z[k][excitation]=circuit_values(names[k])[1]
        own=Z[excitation][excitation]
        check('Потери не ниже DC-границы медных витков '+row['model'],max(0,dc/own.real-1),.002)
        row.update({'matrix_column_complex_ohm':[pair(Z[k][excitation]) for k in range(2)],'winding_R_ohm':own.real,'dc_winding_R_ohm':dc,'primary':abs(own)})
        records.append(row);finish()
    check('Комплексная взаимность катушек '+str(f)+' сетка '+str(level+1),abs(Z[0][1]-Z[1][0])/max(abs(Z[0][1]),abs(Z[1][0])),.003)
    check('Зеркальная симметрия собственных импедансов '+str(f)+' сетка '+str(level+1),abs(Z[0][0]/Z[1][1]-1),.005)
    zc=1/(1j*2*math.pi*f*C)
    a=Z[0][0]+RLEAD+zc;d=Z[1][1]+RLEAD+RLOAD+zc;det=a*d-Z[0][1]*Z[1][0]
    currents=[VOLTAGE*d/det,-VOLTAGE*Z[1][0]/det]
    names=build(f,level,currents)
    row=solve(level,f,'-loaded');row['shape']='refinement-loaded' if refinement else 'loaded';row['excitation']=2
    energy_identity(row,names,f)
    values=[circuit_values(name) for name in names]
    v=[item[1] for item in values]
    residual1=v[0]+(RLEAD+zc)*currents[0]-VOLTAGE
    residual2=v[1]+(RLEAD+RLOAD+zc)*currents[1]
    source=VOLTAGE*currents[0].conjugate()/2
    load=RLOAD*abs(currents[1])**2/2
    leads=RLEAD*sum(abs(ii)**2 for ii in currents)/2
    field=sum(v[k]*currents[k].conjugate()/2 for k in (0,1))
    capacitors=zc*sum(abs(ii)**2 for ii in currents)/2
    total=field+leads+load+capacitors
    eta=load/source.real
    check('Два уравнения Кирхгофа с компенсацией '+row['model'],max(abs(residual1),abs(residual2))/VOLTAGE,.0002)
    check('Комплексный баланс S источника и всех элементов '+row['model'],abs(source-total)/max(abs(source),1e-20),.0002)
    check('Активная мощность — медь, выводы и нагрузка '+row['model'],abs(source.real-row['ohmic_loss_W']-leads-load)/max(source.real,1e-20),.015)
    check('Пассивный КПД между нулём и единицей '+row['model'],max(0,-eta,eta-1),0)
    row.update({'compensation_C_F':C,'frequency_Hz':f,'primary_current_complex_A':pair(currents[0]),'secondary_current_complex_A':pair(currents[1]),'primary_current_amplitude_A':abs(currents[0]),'secondary_current_amplitude_A':abs(currents[1]),'source_active_power_W':source.real,'source_reactive_power_var':source.imag,'field_reactive_power_var':field.imag,'capacitor_reactive_power_var':capacitors.imag,'lead_losses_W':leads,'load_power_W':load,'efficiency':eta,'KVL_relative':max(abs(residual1),abs(residual2))/VOLTAGE,'complex_power_relative':abs(source-total)/max(abs(source),1e-20),'input_impedance_real_ohm':(VOLTAGE/currents[0]).real,'input_impedance_imag_ohm':(VOLTAGE/currents[0]).imag,'primary':load})
    records.append(row);finish();cache[(f,level)]=row;return row

for f in (60000,80000,100000,120000,140000):
    for level in range(3):measure(f,level)
history=[]
for spacing in (20000,10000,5000,2500,1250):
    finest=[r for (f,level),r in cache.items() if level==2]
    best=max(finest,key=lambda row:row['load_power_W'])
    if spacing<20000:
        for f in (best['frequency_Hz']-spacing,best['frequency_Hz']+spacing):measure(f,2,True)
        best=max((r for (f,level),r in cache.items() if level==2),key=lambda row:row['load_power_W'])
    history.append({'frequency_step_Hz':spacing,'peak_frequency_Hz':best['frequency_Hz'],'peak_load_power_W':best['load_power_W'],'peak_efficiency':best['efficiency']})
check('Уточнение частотного максимума мощности 2500→1250Гц',abs(history[-1]['peak_load_power_W']/history[-2]['peak_load_power_W']-1),.015)
check('Максимум лежит внутри проверенного частотного диапазона',max(0,60000-history[-1]['peak_frequency_Hz'],history[-1]['peak_frequency_Hz']-140000),0)
check('Есть разрешимая резонансная селективность',max(0,1.2-history[-1]['peak_load_power_W']/min(cache[(60000,2)]['load_power_W'],cache[(140000,2)]['load_power_W'])),0)
records.append({'measurement':'frequency_refinement','compensation_C_F':C,'frequency_refinement_history':history})
groups={}
for row in records:
    if 'level' not in row or row.get('shape','').startswith('refinement'):continue
    key=str(row['variant'])+row['shape']+str(row['excitation'])
    groups.setdefault(key,[]).append(row)
for key,rows in groups.items():
    rows.sort(key=lambda row:row['level'])
    check('Сходимость по трём сеткам '+key,abs(rows[-1]['primary']/rows[-2]['primary']-1),.02)
result={'runtime':'FEMM 5; embedded MicroPython','case_id':'resonant-wireless-power','records':records,'checks':checks,'passed':all(c['passed'] for c in checks),'scope':'Самостоятельная гармоническая осесимметричная модель двух катушек с последовательной C-компенсацией и пассивной нагрузкой; численные проверки относятся к указанным параметрам.'}
with open(directory+'/result.json','w') as stream:json.dump(result,stream)
print(json.dumps(result));assert result['passed']