Индукционный каротаж: передатчик, приёмник и слоистая проводимость

"""Собственная учебная лабораторная FEMM 5: встроенный Python.
Команды FEMM доступны глобально, внешний пакет femm не используется.
Каждая сетка сохраняется в отдельный FEM, закрывается, открывается повторно
и решается с русскими именами свойств. Новый запуск создаёт calculated-NN.
Опубликованные модели и квитанции не изменяются. Размеры: мм; результаты: СИ.
"""
import os
import sys
import math
import json
CONFIG = {'N': 100, 'I': 1, 'f': 100000, 'variants': [0, 1, 2, 3, 4]}
MODE = 'logging'
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=='lvdt':
 for variant in CONFIG['variants']:
  for excitation in (1,2):
   for level in range(3):
    start(True,0,250,level);material('Медь идеальных обмоток');material('Линейный сердечник LVDT',CONFIG['mu'])
    rect(0,variant-15,3,variant+15);label(1,variant,'Линейный сердечник LVDT',.6/(2**level),5)
    names=[]
    specs=[(0,200,1 if excitation==1 else 0),(20,100,1 if excitation==2 else 0),(-20,100,0)]
    for num,(z,turns,current) in enumerate(specs,1):
        name=('Первичная обмотка','Верхняя вторичная обмотка','Нижняя вторичная обмотка')[num-1]
        mi_addcircprop(name,current,1);names.append(name)
        rect(6,z-4,10,z+4);label(8,z,'Медь идеальных обмоток',.6/(2**level),10+num,name,turns)
    rect(0,-60,30,60);label(20,0,'Воздух',1.2/(2**level),3);label(150,0,'Воздух',15/(2**level),1)
    row=solve(level,variant,'-exc%d'%excitation);row['excitation']=excitation
    energy_identity(row,names,0)
    flux=[circuit_values(n)[2].real for n in names]
    if excitation==1:
        row.update({'primary_L_H':flux[0],'upper_M_H':flux[1],'lower_M_H':flux[2],'differential_M_H':flux[1]-flux[2],'differential_voltage_amplitude_V':2*math.pi*CONFIG['f']*abs(flux[1]-flux[2]),'differential_voltage_signed_V':2*math.pi*CONFIG['f']*(flux[1]-flux[2]),'primary':flux[0]})
    else:row.update({'reverse_upper_M_H':flux[0],'upper_self_L_H':flux[1],'primary':flux[1]})
    records.append(row);finish()
 for displacement in CONFIG['variants']:
    a=[r for r in records if r['variant']==displacement and r['level']==3 and r['excitation']==1][0]
    b=[r for r in records if r['variant']==displacement and r['level']==3 and r['excitation']==2][0]
    check('Взаимность первичная — верхняя при смещении '+str(displacement),abs(a['upper_M_H']/b['reverse_upper_M_H']-1),.002)
    if displacement==0:check('Дифференциальный нуль в центре',abs(a['differential_M_H'])/max(abs(a['upper_M_H']),abs(a['lower_M_H'])),.002)
    else:check('Знак дифференциального отклика при смещении '+str(displacement),max(0,-a['differential_M_H']*displacement),0)
 for displacement in (5,10):
    positive=[r for r in records if r['variant']==displacement and r['level']==3 and r['excitation']==1][0]
    negative=[r for r in records if r['variant']==-displacement and r['level']==3 and r['excitation']==1][0]
    check('Нечётность калибровочной характеристики ±'+str(displacement),abs(positive['differential_M_H']+negative['differential_M_H'])/max(abs(positive['differential_M_H']),abs(negative['differential_M_H'])),.01)

elif MODE=='logging':
 def filament_mutual(radius,spacing,turns):
    count=720;step=2*math.pi/count;total=0
    for k in range(count+1):
        angle=k*step
        value=math.cos(angle)/math.sqrt(2*radius*radius+spacing*spacing-2*radius*radius*math.cos(angle))
        total+=value*(1 if k in (0,count) else 4 if k%2 else 2)
    return MU0*turns*turns*radius*radius/2*total*step/3
 controls=((0,0),(.1,.1),(1,1),(10,10),(.1,10))
 for variant in CONFIG['variants']:
  for excitation in (1,2):
   for level in range(3):
    # Прямоугольная область: r=0..20м,z=±20м; ближние области сгущаются отдельно.
    newdocument(0);mi_probdef(CONFIG['f'],'millimeters','axi',1e-10,1000,30,0)
    material('Воздух скважины');material('Медь идеальных катушек')
    lo,hi=controls[variant]
    material('Нижний пласт',1,lo);material('Верхний пласт',1,hi)
    mi_addboundprop('Удалённая граница A=0',0,0,0,0,0,0,0,0,0)
    polygon([(0,-20000),(75,-20000),(1500,-20000),(20000,-20000),(20000,0),(20000,20000),(1500,20000),(75,20000),(0,20000)],'Удалённая граница A=0',1000/(2**level))
    # Полный цилиндр скважины и два горизонтальных пласта.
    for y in (-1000,0,1000):mi_addnode(75,y)
    for a,b in ((-20000,-1000),(-1000,0),(0,1000),(1000,20000)):mi_addsegment(75,a,75,b)
    mi_addsegment(75,0,20000,0)
    # Ближняя зона горной среды и скважины не изменяет физических материалов.
    for y in (-1000,1000):
        mi_addnode(1500,y);mi_addsegment(75,y,1500,y)
    mi_addnode(1500,0);mi_addsegment(1500,-1000,1500,0);mi_addsegment(1500,0,1500,1000)
    for y in (-750,750):
        mi_addnode(0,y);mi_addnode(75,y);mi_addsegment(0,y,75,y)
    label(40,0,'Воздух скважины',8/(2**level),3)
    label(40,-10000,'Воздух скважины',1000/(2**level),1);label(40,10000,'Воздух скважины',1000/(2**level),1)
    label(500,-500,'Нижний пласт',60/(2**level),5);label(500,500,'Верхний пласт',60/(2**level),6)
    label(5000,-5000,'Нижний пласт',600/(2**level),5);label(5000,5000,'Верхний пласт',600/(2**level),6)
    names=[]
    for num,z in enumerate((-250,250),1):
        name='Передающая катушка' if num==1 else 'Приёмная катушка'
        current=1 if num==excitation else 0
        mi_addcircprop(name,current,1);names.append(name)
        rect(48,z-5,52,z+5);label(50,z,'Медь идеальных катушек',1/(2**level),10+num,name,CONFIG['N'])
    row=solve(level,variant,'-exc%d'%excitation);row['excitation']=excitation
    energy_identity(row,names,CONFIG['f'])
    ii,vv,mm=circuit_values(names[1 if excitation==1 else 0]);transfer=mm
    row.update({'sigma_lower_S_m':lo,'sigma_upper_S_m':hi,'mutual_real_H':transfer.real,'mutual_imag_H':transfer.imag,'mutual_abs_H':abs(transfer),'transfer_phase_deg':math.atan2(transfer.imag,transfer.real)*180/math.pi,'receiver_voltage_amplitude_V':2*math.pi*CONFIG['f']*abs(transfer),'primary':abs(transfer)})
    if variant==0:
        exact=filament_mutual(.05,.5,CONFIG['N']);row['Neumann_air_mutual_H']=exact
        check('Независимый интеграл Неймана в воздухе '+row['model'],abs(transfer.real/exact-1),.025)
        check('Воздушный контроль без фазового сдвига '+row['model'],abs(transfer.imag)/abs(transfer),1e-8)
    records.append(row);finish()
 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]
    ma=complex(a['mutual_real_H'],a['mutual_imag_H']);mb=complex(b['mutual_real_H'],b['mutual_imag_H'])
    check('Комплексная взаимность Tx/Rx в среде '+str(variant),abs(ma-mb)/abs(ma),.003)
 air=[r for r in records if r['variant']==0 and r['level']==3 and r['excitation']==1][0]
 homogeneous=[r for r in records if r['variant']==2 and r['level']==3 and r['excitation']==1][0]
 layered=[r for r in records if r['variant']==4 and r['level']==3 and r['excitation']==1][0]
 high=[r for r in records if r['variant']==3 and r['level']==3 and r['excitation']==1][0]
 signal=abs(complex(high['mutual_real_H']-air['mutual_real_H'],high['mutual_imag_H']-air['mutual_imag_H']))/air['mutual_abs_H']
 layer_signal=abs(complex(layered['mutual_real_H']-homogeneous['mutual_real_H'],layered['mutual_imag_H']-homogeneous['mutual_imag_H']))/homogeneous['mutual_abs_H']
 check('Разрешимый отклик на проводимость среды',max(0,.01-signal),0)
 check('Разрешимый отклик на слоистость',max(0,.005-layer_signal),0)

elif MODE=='meissner':
 for variant in CONFIG['variants']:
  for level in range(3):
    radius=variant;outer=20*radius
    start(True,0,outer,level,CONFIG['B'])
    mi_addboundprop('Идеальная граница Мейснера',0,0,0,0,0,0,0,0,0)
    halfcircle(radius,2/(2**level),'Идеальная граница Мейснера')
    label(radius*.5,0,'<No Mesh>',radius/10,0)
    halfcircle(3*radius,3/(2**level))
    label(2*radius,0,'Воздух',radius/10/(2**level),3);label(10*radius,0,'Воздух',outer/20/(2**level),1)
    row=solve(level,variant)
    probes=[]
    for rp,zp in ((1.5,0),(2,0),(0,1.5),(0,2),(.9,.9)):
        rr=math.sqrt(rp*rp+zp*zp);alpha=1/rr**3
        br=-1.5*CONFIG['B']*alpha*rp*zp/(rr*rr)
        bz=CONFIG['B']*((1-alpha)*zp*zp+(1+alpha/2)*rp*rp)/(rr*rr)
        measured_r,measured_z=magnetic_b(rp*radius,zp*radius)
        measured_r=real(measured_r);measured_z=real(measured_z)
        err=math.sqrt((measured_r-br)**2+(measured_z-bz)**2)/CONFIG['B']
        check('Аналитическое поле идеальной диамагнитной сферы '+str(rp)+','+str(zp)+' '+row['model'],err,.035/(2**level))
        probes.append([rp,zp,measured_r,measured_z,br,bz])
    row.update({'radius_mm':radius,'equator_1p5a_Bz_T':probes[0][3],'pole_1p5a_Bz_T':probes[2][3],'field_probes_r_over_a_z_over_a_T':probes,'energy_per_radius_cubed_J_m3':row['energy_J']/(radius*.001)**3,'primary':probes[0][3]})
    records.append(row);finish()
 finest=[r for r in records if r['level']==3]
 for row in finest[1:]:check('Геометрическое подобие энергии a='+str(row['radius_mm']),abs(row['energy_per_radius_cubed_J_m3']/finest[0]['energy_per_radius_cubed_J_m3']-1),.005)
# Сходимость по трём сеткам проверяется отдельно для каждого физического режима.
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 ('lvdt','logging','meissner') else 0.012)
result={'runtime':'FEMM 5; embedded MicroPython','case_id':'electromagnetic-well-logging','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']