"""Силы в сохранённом полном периоде: встроенный MicroPython FEMM 5. Сначала выполните этап15 (или all) из reproduce.py. Задайте TG_RESULT_FILE абсолютным именем его полного result.json. Все144 FEM/ANS должны лежать рядом. FORCE_KIND='both'|'healthy'|'eccentric'; FORCE_SLOT=6 — номер паза от1. FORCE_POSITIONS=None — полный период; [0,18,36,54] — публичные контрольные точки72. TG_FORCE_OUTPUT_ROOT — существующий каталог для нового forces-calculated-NN. Входы открываются только для чтения: mi_analyze/mi_saveas здесь отсутствуют. Сохранённые публичные контрольные точки не заменяют полный расчёт этапа15. """ import os import sys import math import json import hashlib import binascii saved=globals().get('__file__','') or (sys.argv[0] if sys.argv else '') SCRIPT=globals().get('TG_FORCE_SCRIPT_FILE',os.path.abspath(saved)) BASE=os.path.dirname(SCRIPT) INPUT=globals().get('TG_RESULT_FILE','') if not INPUT: raise ValueError('Задайте TG_RESULT_FILE: полный result.json после STAGE15 или all') SOURCE=os.path.dirname(os.path.abspath(INPUT)) DEST=globals().get('TG_FORCE_OUTPUT_ROOT',SOURCE) KIND=globals().get('FORCE_KIND','both') SLOT=int(globals().get('FORCE_SLOT',6)) POSITIONS=globals().get('FORCE_POSITIONS',None) ENGINE=globals().get('TG_ENGINE_FILE','C:/Program Files/femm-qt/bin/femm_headless.exe') ENGINE_VERSION=globals().get('TG_ENGINE_VERSION','FEMM headless engine v5.0 06Oct2026') def sha(path): h=hashlib.sha256() with open(path,'rb') as f: while True: data=f.read(65536) if not data: break h.update(data) return binascii.hexlify(h.digest()).decode() def real(value): return float(complex(value).real) with open(INPUT) as f: original=json.load(f) if '15' not in original.get('results',{}): raise ValueError('Во входном result.json отсутствует завершённый этап15') if KIND not in ['both','healthy','eccentric']: raise ValueError('FORCE_KIND: both, healthy или eccentric') if POSITIONS is not None and not POSITIONS: raise ValueError('FORCE_POSITIONS: непустой список индексов либо None') cfg=original['parameters'] if original.get('geometry_revision')!='open-slot-wedges-v2' or original.get('geometry_contract',{}).get('revision')!='open-slot-wedges-v2': raise ValueError('Нужен результат новой открытой геометрии open-slot-wedges-v2') if not 1<=SLOT<=cfg['stator_slots']: raise ValueError('FORCE_SLOT вне диапазона пазов статора') kinds=['healthy','eccentric'] if KIND=='both' else [KIND] angle=(SLOT-.5)*360/cfg['stator_slots'] radius=sum(cfg['stator_slot_radii_mm'])/2 x,y=radius*math.cos(angle*math.pi/180),radius*math.sin(angle*math.pi/180) phase,sign=[(0,1),(2,-1),(1,1),(0,-1),(2,1),(1,-1)][int(((angle-60)%360)/60)] circuit=['A — фаза статора','B — фаза статора','C — фаза статора'][phase] turns=sign*cfg['slot_turns_stator'] profile=original['geometry_contract']['profiles']['stator'] # Медь нового открытого паза — прямоугольник: S=bs·(xg−xe), мм²→м². # Площадь совпадает со старой лишь по условию сравнения; топология паза изменилась. area=profile['bs_mm']*abs(profile['xg_mm']-profile['xe_mm'])*1e-6 if abs(area-profile['copper_area_mm2']*1e-6)/area>1e-10: raise ValueError('Прямоугольник меди не согласован с опубликованным контрактом') volume=area*cfg['depth_mm']*.001 index=1 while os.path.exists(DEST+'/forces-calculated-%02d'%index): index+=1 OUT=DEST+'/forces-calculated-%02d'%index os.mkdir(OUT) records=[] checks=[] inputs=[] for kind in kinds: wave=original['results']['15'][kind] count=len(wave['samples']) if POSITIONS is not None and any(int(k)!=k or k<0 or k>=count for k in POSITIONS): raise ValueError('FORCE_POSITIONS вне диапазона сохранённых положений') for k,source in enumerate(wave['samples']): if POSITIONS is not None and k not in POSITIONS: continue model=SOURCE+'/'+source['model'] solution=SOURCE+'/'+source['solution'] if not os.path.exists(model) or not os.path.exists(solution): raise ValueError('Нужен полный сырой результат этапа15; отсутствует '+source['solution']) model_sha,solution_sha=sha(model),sha(solution) opendocument(model) mi_loadsolution() current=real(mo_getcircuitproperties(circuit)[0]) mo_clearblock() mo_selectblock(x,y) # Однородный DC-ток проводника: Jz=N·Iф/S. Пассивный ток FEMM уже имеет знак. jz=turns*current/area native_area=real(mo_blockintegral(5)) native_volume=real(mo_blockintegral(10)) native_current=real(mo_blockintegral(7)) bx,by=real(mo_blockintegral(8)),real(mo_blockintegral(9)) fx,fy=real(mo_blockintegral(11)),real(mo_blockintegral(12)) # Jz·ez×B=(−Jz·By,+Jz·Bx); интегралы8/9 уже включают активную длину. reconstructed=[-jz*by,jz*bx] identity=math.sqrt((fx-reconstructed[0])**2+(fy-reconstructed[1])**2)/max(math.sqrt(fx*fx+fy*fy),1e-12) area_error=abs(native_area-area)/area volume_error=abs(native_volume-volume)/volume current_error=abs(native_current-turns*current)/max(abs(turns*current),1.0) mo_clearblock() mo_groupselectblock(2) total=[real(mo_blockintegral(18)),real(mo_blockintegral(19))] conductors=[real(mo_blockintegral(11)),real(mo_blockintegral(12))] iron=[total[j]-conductors[j] for j in range(2)] total_torque=real(mo_blockintegral(22)) conductor_torque=real(mo_blockintegral(15)) # Суммарная сила тела=сила на токах+сила на намагниченном железе. # Вал и клинья немагнитны (μr=1) и не несут J; разность — сила железа. reread_error=max(abs(total[j]-source['rotor_force_N'][j]) for j in range(2))/max(1.0,math.sqrt(sum(v*v for v in total))) mo_clearblock() mo_close() mi_close() if sha(model)!=model_sha or sha(solution)!=solution_sha: raise RuntimeError('Постобработка изменила входной файл') row={'kind':kind,'position_index':k,'mechanical_angle_deg':k*360/count, 'time_s':k/(cfg['frequency_Hz']*count),'source_model':source['model'], 'source_solution':source['solution'],'source_model_sha256':model_sha,'source_solution_sha256':solution_sha, 'slot_number':SLOT,'slot_center_mm':[x,y],'circuit':circuit,'signed_turns':turns, 'phase_current_A':current,'uniform_Jz_A_m2':jz,'slot_area_m2':area,'slot_volume_m3':volume, 'native_area_m2':native_area,'native_volume_m3':native_volume,'native_total_current_Aturn':native_current, 'integral_B_T_m3':[bx,by],'native_lorentz_N':[fx,fy],'reconstructed_force_N':reconstructed, 'relative_identity_error':identity,'area_relative_error':area_error, 'volume_relative_error':volume_error,'current_relative_error':current_error, 'rotor_total_weighted_N':total,'rotor_conductor_lorentz_N':conductors, 'rotor_ferromagnetic_difference_N':iron,'rotor_total_torque_origin_Nm':total_torque, 'rotor_conductor_torque_origin_Nm':conductor_torque, 'rotor_ferromagnetic_torque_difference_origin_Nm':total_torque-conductor_torque, 'reread_force_relative_error':reread_error,'inputs_unchanged':True} records.append(row) print('FORCE',kind,k,fx,fy,identity) checks=[{'name':'Тождество Лоренца J×B','passed':all(r['relative_identity_error']<1e-6 for r in records),'limit':1e-6, 'max_relative_error':max(r['relative_identity_error'] for r in records)}, {'name':'Площадь и объём выбранного паза','passed':all(r['area_relative_error']<1e-8 and r['volume_relative_error']<1e-8 for r in records),'limit':1e-8}, {'name':'Ампер-витки N·I и интегралJ','passed':all(r['current_relative_error']<1e-8 for r in records),'limit':1e-8}, {'name':'Повторное чтение силы ротора','passed':all(r['reread_force_relative_error']<1e-10 for r in records),'limit':1e-10}, {'name':'Все входные FEM/ANS неизменны','passed':all(r['inputs_unchanged'] for r in records)}] result={'schema':'femm5-saved-conductor-forces-v2','geometry_revision':original['geometry_revision'], 'copper_area_formula':'S_m2 = bs_mm * abs(xg_mm - xe_mm) * 1e-6', 'runtime':sys.version,'implementation':str(sys.implementation), 'engine_version':ENGINE_VERSION,'engine_file_sha256':sha(ENGINE),'script_sha256':sha(SCRIPT), 'input_result_sha256':sha(INPUT),'input_result_file':os.path.abspath(INPUT), 'force_kind':KIND,'slot_number':SLOT,'records':records,'checks':checks, 'position_selection':'all' if POSITIONS is None else POSITIONS, 'full_period_processed':POSITIONS is None, 'passed':all(c['passed'] for c in checks), 'scope':'Постобработка точных сохранённых 2D-решений. Лоренц для одного однородного проводника; сила железа ротора как разность агрегированной взвешенной силы тела и Лоренца всех его токов. Это проверка внутренней согласованности интегралов, не независимая полевая/физическая проверка или сходимость сил; торцевые и механические нагрузки не моделируются.'} with open(OUT+'/conductor-forces.json','w') as f: json.dump(result,f) with open(OUT+'/conductor-forces.csv','w') as f: f.write('kind,angle_deg,time_s,phase_current_A,Fx_lorentz_N,Fy_lorentz_N,Fx_JcrossB_N,Fy_JcrossB_N,identity_relative_error,Fx_rotor_total_N,Fy_rotor_total_N,Fx_rotor_conductors_N,Fy_rotor_conductors_N,Fx_rotor_iron_difference_N,Fy_rotor_iron_difference_N\n') for r in records: values=[r['kind'],r['mechanical_angle_deg'],r['time_s'],r['phase_current_A']]+r['native_lorentz_N']+r['reconstructed_force_N']+[r['relative_identity_error']]+r['rotor_total_weighted_N']+r['rotor_conductor_lorentz_N']+r['rotor_ferromagnetic_difference_N'] f.write(','.join(str(v) for v in values)+'\n') print('FORCE_RESULT',OUT,result['passed'],len(records)) if not result['passed']: raise RuntimeError('Проверки сил не пройдены; изучите новый conductor-forces.json')