# Execute with FEMM5 native embeddedMicroPython from this extracted folder. import os def _exists(p): try:os.stat(p);return True except OSError:return False if _exists('files/reproduce.py'):os.chdir('files') elif '__file__' in globals() and '/' in __file__.replace('\\','/'):os.chdir(__file__.replace('\\','/').rsplit('/',1)[0]) import json,os,sys,time,math,cmath import numpy as np CONFIG={'output': 'native-circuits-native-result.json', 'directories': {'practice-ii04-1': '.', 'practice-ii03-3': '.'}, 'models': {'./practice-ii04-1-source-model.fem': {'problem': 'planar', 'depth': 40.0, 'angle': 25.0, 'circuits': ['Первичная фаза А', 'Первичная фаза В', 'Первичная фаза С', 'Вторичная фаза А', 'Вторичная фаза В', 'Вторичная фаза С']}, './secondary-marking-reversed-source-model.fem': {'problem': 'planar', 'depth': 40.0, 'angle': 25.0, 'circuits': ['Первичная фаза А', 'Первичная фаза В', 'Первичная фаза С', 'Вторичная фаза А', 'Вторичная фаза В', 'Вторичная фаза С']}, './practice-ii03-3-source-model.fem': {'problem': 'planar', 'depth': 40.0, 'angle': 25.0, 'circuits': ['Первичная обмотка', 'Вторичная обмотка']}, './strict-nonlinear-source-model.fem': {'problem': 'planar', 'depth': 40.0, 'angle': 25.0, 'circuits': ['Первичная обмотка', 'Вторичная обмотка']}}} DATA={'runtime':sys.version,'implementation':str(sys.implementation),'experiments':{},'field_solutions':0} assert sys.implementation.name=='micropython' def save(): with open(CONFIG['output'],'w') as f:json.dump(DATA,f) def serial(v): if isinstance(v,complex):return [v.real,v.imag] if isinstance(v,(tuple,list)):return [serial(x) for x in v] if isinstance(v,dict):return {k:serial(x) for k,x in v.items()} return v def close(): for n in ('mo_close','mi_close'): try:globals()[n]() except Exception:pass def start(name): close();d=CONFIG['directories'][name];os.chdir(d);return d def model(name):return './'+name+'-source-model.fem' def solve(path,tag,currents=None,frequency=None,rotate=0,group=2,mesh=None,material=None): close();opendocument(path);mi_saveas(os.getcwd()+'/'+tag+'.fem') if frequency is not None: info=CONFIG['models'][path] mi_probdef(frequency,'millimeters',info['problem'],1e-8,info['depth'],info['angle']) if rotate: mi_selectgroup(group);mi_moverotate(0,0,rotate);mi_clearselected() if material: x,y,material_name=material;mi_selectlabel(x,y);mi_setblockprop(material_name,0,1,'',0,0,0);mi_clearselected() if mesh: for x,y,mat,sz,circ,mag,g,turns in mesh: mi_selectlabel(x,y);mi_setblockprop(mat,0,sz,circ,mag,g,turns);mi_clearselected() if currents: for name,value in currents:mi_modifycircprop(name,1,value) mi_analyze(1);mi_loadsolution();DATA['field_solutions']+=1 info={'file':os.getcwd()+'/'+tag+'.fem','nodes':mo_numnodes(),'elements':mo_numelements(),'circuits':{}} for name in CONFIG['models'][path]['circuits']: info['circuits'][name]=list(mo_getcircuitproperties(name)) mo_clearblock();mo_groupselectblock();info['energy_J']=mo_blockintegral(2);info['coenergy_J']=mo_blockintegral(17) info['loss_magnetic_W']=mo_blockintegral(3);info['loss_resistive_W']=mo_blockintegral(4);info['loss_total_W']=mo_blockintegral(6) mo_clearblock() try: mo_groupselectblock(group);info['torque_Nm']=mo_blockintegral(22);info['force_x_N']=mo_blockintegral(18);info['force_y_N']=mo_blockintegral(19);mo_clearblock() except Exception as error:info['force_error']=str(error) return info def real(x):return complex(x).real def add(name,value):DATA['experiments'][name]=serial(value);save() def z_matrix(path,prefix,names): columns=[];evidence=[] for j in range(len(names)): row=solve(path,prefix+'-basis-'+str(j),[(n,1 if j==k else 0) for k,n in enumerate(names)]) columns.append([complex(row['circuits'][n][1]) for n in names]);evidence.append(row) return [[columns[j][i] for j in range(len(names))] for i in range(len(names))],evidence def linear_load(path,prefix,names,Z,rloads,rprimary,source,neutral=False): n=len(names);m=n+1 if neutral else n matrix=[[0j for j in range(m)] for i in range(m)] rhs=[0j]*m for i in range(n): for j in range(n):matrix[i][j]=Z[i][j] if i=n//2:r+=complex(answer[n]) residual.append(r) row['current_A']=currents;row['voltage_V']=voltages;row['residual_V']=residual row['source_power_W']=sum(real(source[i]*currents[i].conjugate())/2 for i in range(n//2)) row['external_primary_copper_W']=sum(rprimary*abs(currents[i])**2/2 for i in range(n//2)) row['load_power_W']=sum(rloads[i-n//2]*abs(currents[i])**2/2 for i in range(n//2,n)) row['neutral_current_A']=-sum(currents[n//2:]);row['neutral_potential_V']=complex(answer[n]) if neutral else 0 row['power_balance_W']=row['source_power_W']-row['external_primary_copper_W']-row['load_power_W']-real(row['loss_total_W']) return row def restore(z):return complex(z[0],z[1]) if isinstance(z,list) else complex(z) def load_previous(): with open('native-experiments-02.json') as f:return json.load(f)['experiments'] try: old=load_previous() # Correct Y-delta connection: each secondary winding is a delta branch, # D maps branch currents to terminal line currents. Each terminal load # is R in star; port matrix added to the FEMM winding voltages is R D^T D. name='practice-ii04-1';start(name);path=model(name);names=CONFIG['models'][path]['circuits'];Z=[[restore(z) for z in row] for row in old['three_phase']['Z_ohm']] D=[[1,0,-1],[-1,1,0],[0,-1,1]];R=10.;matrix=[[Z[i][j] for j in range(6)] for i in range(6)] for i in range(3):matrix[i][i]+=.8 for i in range(3): for j in range(3):matrix[i+3][j+3]+=R*sum(D[k][i]*D[k][j] for k in range(3)) source=[400*math.sqrt(2)/math.sqrt(3)*cmath.exp(-2j*math.pi*k/3) for k in range(3)] solution=np.linalg.solve(np.array(matrix),np.array(source+[0j]*3));I=[complex(v) for v in solution] r=solve(path,'verified-y-delta',list(zip(names,I)));v=[complex(r['circuits'][n][1]) for n in names] line_currents=[sum(D[k][j]*I[j+3] for j in range(3)) for k in range(3)] branch_residual=[v[i+3]+R*sum(D[k][i]*line_currents[k] for k in range(3)) for i in range(3)] r['delta_incidence_D']=D;r['winding_currents_amplitude_A']=I;r['secondary_line_currents_amplitude_A']=line_currents;r['secondary_winding_voltage_amplitude_V']=v[3:];r['branch_kvl_residual_V']=branch_residual;r['primary_kvl_residual_V']=[v[i]+.8*I[i]-source[i] for i in range(3)];r['secondary_line_voltage_amplitude_V']=v[3:];r['primary_line_voltage_amplitude_V']=[source[0]-source[1],source[1]-source[2],source[2]-source[0]];r['source_power_W']=sum(real(source[k]*I[k].conjugate())/2 for k in range(3));r['load_power_W']=R*sum(abs(x)**2 for x in line_currents)/2;r['primary_copper_W']=.8*sum(abs(x)**2 for x in I[:3])/2;r['balance_W']=r['source_power_W']-r['load_power_W']-r['primary_copper_W']-real(r['loss_total_W']) add('star_delta_connection',r) # Nonsingular source/load operating points with an actual nonlinear steel # law. Every evaluation of residual is a freshly solved FEMM field problem. name='practice-ii03-3';start(name);path=model(name);names=CONFIG['models'][path]['circuits'];calls=0 def residual(x,U,hlabel): global calls currents=[complex(x[0],x[1]),complex(x[2],x[3])] item=solve(path,'newton-transient-current',list(zip(names,currents)),material=(0,40,'Нелинейная сталь')) calls+=1 item.pop('file');item['replay_source']=path;item['requested_currents_A']=currents u=[complex(item['circuits'][n][1]) for n in names] z=[u[0]+.8*currents[0]-U,u[1]+10*currents[1]] return [z[0].real,z[0].imag,z[1].real,z[1].imag],item def norm(x):return max(abs(v) for v in x) def newton(U,initial,h): global calls initial_calls=calls;x=list(initial);log=[];success=False for iteration in range(20): r,field=residual(x,U,'base');scale=norm(r)/max(abs(U),1.) entry={'iteration':iteration,'current_components_A':list(x),'residual_components_V':r,'relative_residual':scale,'field_solve_count':calls-initial_calls} log.append(entry) if scale<1e-5:success=True;break columns=[] for j in range(4): xp=list(x);xm=list(x);xp[j]+=h;xm[j]-=h;rp,_=residual(xp,U,'plus');rm,_=residual(xm,U,'minus');columns.append([(rp[i]-rm[i])/(2*h) for i in range(4)]) J=[[columns[j][i] for j in range(4)] for i in range(4)];dx=np.linalg.solve(np.array(J),np.array([-v for v in r]));alpha=1. while alpha>=1/64: candidate=[x[j]+alpha*float(dx[j]) for j in range(4)] if max(abs(v) for v in candidate)>100:alpha/=2;continue test,_=residual(candidate,U,'trial') if norm(test)