"""Учебный турбогенератор: 16 воспроизводимых этапов, FEMM 5 / MicroPython. Запускать внутри FEMM 5. Команды mi_*/mo_* доступны глобально. STAGE = 'all' либо '01'..'16'; CASE = 'standard' (или 'smoke'). TG_ROOT можно задать абсолютным каталогом этого файла перед exec(). Каждый запуск создаёт calculated-NN. Исходники и старые результаты сохраняются. Геометрия и B(H) авторские, условные; это не копия ТГ 340 МВт из статей. Открытые пазы со шлицем и подклиновым расширением; отдельные условные μr=1 клинья. Профиль статора следует опубликованному Lua (Милых, Полякова, 2015, с.11), размеры и граница воздух/клин авторские; профиль ротора — зеркальная адаптация. Цилиндрический ротор, 2D активная длина, без демпфера, торцов, вихревых токов, механики и физической аттестации реальной машины. """ import os import sys import math import json STAGE = str(globals().get('STAGE', 'all')) CASE = str(globals().get('CASE', 'standard')) saved = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '') BASE = globals().get('TG_ROOT', os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd()) OUTPUT_BASE=globals().get('TG_OUTPUT_ROOT',BASE) index = 1 while os.path.exists(OUTPUT_BASE + '/calculated-%02d' % index): index += 1 OUT = OUTPUT_BASE + '/calculated-%02d' % index os.mkdir(OUT) PI = math.pi FREQ = 50.0 OMEGA = 2 * PI * FREQ MU0 = 4e-7 * PI SQRT2 = math.sqrt(2) PHASES = ['A — фаза статора', 'B — фаза статора', 'C — фаза статора'] FIELD = 'Обмотка возбуждения ротора' BH = [[0, 0], [0.4, 80], [0.8, 160], [1.2, 300], [1.45, 600], [1.6, 1400], [1.75, 4500], [1.9, 14000], [2.1, 50000], [2.3, 150000]] CFG = {'rotor_radius_mm': 70, 'bore_radius_mm': 72, 'stator_radius_mm': 120, 'boundary_radius_mm': 150, 'shaft_radius_mm': 15, 'depth_mm': 250, 'stator_slots': 24, 'rotor_slots': 12, 'slot_turns_stator': 12, 'slot_turns_field': 12, 'stator_slot_radii_mm': [78, 100], 'rotor_slot_radii_mm': [52, 66], 'stator_slot_width_deg': 8, 'rotor_slot_width_deg': 7, 'frequency_Hz': FREQ, 'pole_pairs': 1, 'geometry_revision':'open-slot-wedges-v2', 'stator_neck_height_mm':1, 'rotor_neck_height_mm':1, 'stator_groove_expansion_each_side_mm':1, 'rotor_groove_expansion_each_side_mm':1, 'stator_chamfer_angles_deg':[45,45], 'rotor_chamfer_angles_deg':[45,45], 'yoke_split_radius_mm': 104, 'rotor_slot_angles_deg': [40,60,80,100,120,140,220,240,260,280,300,320], 'phase_resistance_ohm': 0.15, 'target_phase_voltage_rms_V': 80, 'target_current_rms_A': 10, 'target_power_factor': 0.8, 'baseline_excitation_A':30, 'occ_excitation_A':[0,10,20,30,45,60], 'scc_excitation_A':[10,20,30,40], 'load_currents_A':[0,5,10,15,20], 'regulation_currents_A':[0,5,10,15], 'iteration_initial_reactance_ohm':15, 'iteration_initial_excitation_A_per_V':.18, 'eccentricity_mm': 0.4, 'BH_T_A_per_m': BH, 'geometry_origin': 'Авторские размеры; профиль открытого статора по опубликованному Lua 2015, с.11; ротор зеркально адаптирован. Не промышленный прототип.', 'wedge_origin':'Отдельная граница x=xb и материал клина μr=1, σ=0 — учебное допущение. В опубликованном Lua шлиц и подклиновая зона составляют общую область μ0 зазора без отдельной метки клина.', 'BH_origin': 'Условная изотропная безгистерезисная B(H), не паспорт марки стали.'} CFG['BH_yoke_T_A_per_m'] = [[b,h] for b,h in BH] CFG['BH_teeth_T_A_per_m'] = [[b,h*1.1] for b,h in BH] CFG['BH_rotor_T_A_per_m'] = [[b,h*1.2] for b,h in BH] for key,value in globals().get('TG_CONFIG',{}).items(): if key not in CFG: raise ValueError('Неизвестный параметр: '+key) CFG[key]=value FREQ=CFG['frequency_Hz'] OMEGA=2*PI*FREQ if CFG['pole_pairs']!=1 or CFG['stator_slots']%6!=0: raise ValueError('Текущая топология: p=1 и число пазов статора кратно6') if CFG['geometry_revision']!='open-slot-wedges-v2' or CFG['rotor_slots']!=len(CFG['rotor_slot_angles_deg']): raise ValueError('Требуются revision open-slot-wedges-v2 и rotor_slots=len(rotor_slot_angles_deg)') if not (CFG['shaft_radius_mm']0 and expansion>0 and direction*(xc-xb)>0 and direction*(xd-xc)>0 and direction*(xe-xd)>0 and direction*(xg-xe)>0): raise ValueError('Профиль '+kind+': требуется последовательность a,b,c,d,e,g и положительные фаски') wall=[[xa,ya],[xb,ya],[xc,yb],[xd,yb],[xe,ya],[xf,ya],[xg,ya]] mirror=lambda points:[[x,-y] for x,y in reversed(points)] cavity=wall+mirror(wall) wedge=wall[1:5]+mirror(wall[1:5]) copper=[[xe,ya],[xg,ya],[xg,-ya],[xe,-ya]] neck_polygon=[[xa,ya],[xb,ya],[xb,-ya],[xa,-ya]] mouth_half_angle=math.asin(ya/radius) # Хорда aa' не совпадает с окружностью. Для статора вычитаем круговой # сегмент из площади удалённого железа; для ротора этот сегмент прибавляется. cap=.5*radius*radius*(2*mouth_half_angle-math.sin(2*mouth_half_angle)) cap_sign=-1 if stator else 1 neck_area=area(neck_polygon)+cap_sign*cap cavity_area=area(cavity)+cap_sign*cap if neck_area<=0 or abs(cavity_area-area(copper)-area(wedge)-neck_area)>1e-8: raise ValueError('Не замкнулся баланс площадей '+kind) if stator and max(math.sqrt(x*x+y*y) for x,y in cavity)>=CFG['yoke_split_radius_mm']: raise ValueError('Паз статора пересекает спинку') if not stator and min(x for x,y in cavity)<=CFG['shaft_radius_mm']: raise ValueError('Паз ротора пересекает вал') values={'kind':kind,'direction':direction,'surface_radius_mm':radius, 'bs_mm':bs,'bks_mm':2*yb,'b1_mm':expansion,'as1_deg':angles[0],'as2_deg':angles[1], 'xa_mm':xa,'xb_mm':xb,'xc_mm':xc,'xd_mm':xd,'xe_mm':xe,'xf_mm':xf,'xg_mm':xg, 'h_slot_mm':abs(xg-xa),'h_neck_mm':neck,'h_wedge_zone_mm':abs(xe-xa), 'mouth_half_angle_deg':mouth_half_angle*180/PI, 'wall_local_mm':wall,'cavity_chord_polygon_local_mm':cavity, 'neck_chord_polygon_local_mm':neck_polygon,'wedge_polygon_local_mm':wedge, 'copper_polygon_local_mm':copper,'circular_segment_area_mm2':cap, 'circular_segment_sign':cap_sign,'neck_area_mm2':neck_area, 'wedge_area_mm2':area(wedge),'copper_area_mm2':area(copper), 'cavity_removed_iron_area_mm2':cavity_area,'old_closed_copper_area_mm2':old_area, 'copper_area_preserved_relative_error':abs(area(copper)-old_area)/old_area, 'copper_volume_m3':area(copper)*CFG['depth_mm']*1e-9, 'wedge_volume_m3':area(wedge)*CFG['depth_mm']*1e-9, 'probes_local_mm':{'mouth':[radius+direction*.05,0], 'neck':[(radius+xb)/2,0],'wedge':[(xb+xe)/2,0], 'copper':[(xe+xg)/2,0]}} return values PROFILES={'stator':profile('stator'),'rotor':profile('rotor')} GEOMETRY_CONTRACT={'revision':CFG['geometry_revision'],'length_unit':'mm', 'profiles':PROFILES,'stator_slot_centres_deg':[(k+.5)*360/CFG['stator_slots'] for k in range(CFG['stator_slots'])], 'rotor_slot_centres_local_deg':CFG['rotor_slot_angles_deg'],'rotor_centre_mm_formula':'[eccentricity_mm,0]', 'rotation_formula':'[x cos(a)−y sin(a)+ecc, x sin(a)+y cos(a)]', 'surface_arcs':'Только между устьями; дуги через устья отсутствуют. Шлицы связаны с воздушным зазором.', 'source_scope':'Статорный профиль по опубликованному Lua 2015, с.11. Численные размеры, зеркало ротора и граница x=xb между воздухом и клином — самостоятельная учебная адаптация.', 'wedge_materials':{'stator':'Клин статора — немагнитный, условный','rotor':'Клин ротора — немагнитный, условный'}, 'wedge_relative_permeability':1,'wedge_conductivity_MSm':0, 'wedge_scope':'μr=1, σ=0 — условный магнитостатический заменитель; реальные клинья могут быть магнитными и электропроводными.', 'groups':{'статорное железо':1,'весь ротор включая клинья и обмотку':2,'воздушный зазор':3, 'внешний воздух':4,'статорные проводники A B C':[10,11,12],'клинья статора':13}, 'copper_area_scope':'Старая полярная форма заменена прямоугольником при той же площади S; NI/S сохранено. Это контроль сравнения, не сохранение старой топологии.'} GEOMETRY_CONTRACT['legacy_width_parameter_scope']='stator_slot_width_deg и rotor_slot_width_deg задают только площадь прежнего полярного Cu-четырёхугольника для сравнения S. Реальная ширина прямоугольного тела нового паза — bs_mm; ширина устья определяется asin(bs/(2R)).' for kind in ['stator','rotor']: p=PROFILES[kind] centers=GEOMETRY_CONTRACT['stator_slot_centres_deg'] if kind=='stator' else CFG['rotor_slot_angles_deg'] centers=sorted(centers) separation=min((centers[(k+1)%len(centers)]-centers[k])%360 for k in range(len(centers))) if separation<=2*p['mouth_half_angle_deg']: raise ValueError('Устья пересекаются: '+kind) p['minimum_tooth_surface_angle_deg']=separation-2*p['mouth_half_angle_deg'] with open(OUT+'/geometry-contract.json','w') as f: json.dump(GEOMETRY_CONTRACT,f) DRAW_NODES={} DRAW_SEGMENTS={} def node(point): # Единые координаты узла устраняют дубли на границе клин/медь и дуга/устье. key=(round(point[0],8),round(point[1],8)) if key not in DRAW_NODES: mi_addnode(key[0],key[1]) DRAW_NODES[key]=True return key def segment(a,b): a,b=node(a),node(b) key=tuple(sorted([a,b])) if a!=b and key not in DRAW_SEGMENTS: mi_addsegment(a[0],a[1],b[0],b[1]) DRAW_SEGMENTS[key]=True def polyline(points): for k in range(len(points)-1): segment(points[k],points[k+1]) def polygon(points): polyline(points+[points[0]]) def circle(radius, eccentricity=0, boundary=False, arc_deg=2, mouths=None, half_angle=0): # Граница зубцов прерывается у открытых пазов. Сплошная окружность через # устье создала бы ложную железную перемычку; здесь она вообще не рисуется. cuts=[0,90,180,270] mouths=mouths or [] for a in mouths: cuts.extend([(a-half_angle)%360,(a+half_angle)%360]) cuts=sorted(set(round(a%360,10) for a in cuts)) for k in range(len(cuts)): start=cuts[k] stop=cuts[(k+1)%len(cuts)] if k+1= 4 and nonlinear: for material,key in [('Сталь спинки статора — условная','BH_yoke_T_A_per_m'), ('Сталь зубцов статора — условная','BH_teeth_T_A_per_m'), ('Сталь ротора — условная','BH_rotor_T_A_per_m')]: for b,h in CFG[key]: mi_addbhpoint(material,b,h) mi_addboundprop('Внешняя граница A=0', 0, 0, 0, 0, 0, 0, 0, 0, 0) arc_deg = [3, 2, 1][mesh] steel_mesh = [6, 4, 2.5][mesh] slot_mesh = [3, 2, 1.3][mesh] gap_mesh = [1.1, 0.7, 0.45][mesh] circle(outer, boundary=True, arc_deg=arc_deg) stator=CFG['stator_radius_mm'] bore=CFG['bore_radius_mm'] rotor=CFG['rotor_radius_mm'] shaft=CFG['shaft_radius_mm'] split=CFG['yoke_split_radius_mm'] circle(stator, arc_deg=arc_deg) circle(split, arc_deg=arc_deg) stator_centres=GEOMETRY_CONTRACT['stator_slot_centres_deg'] if level>=2 else [] rotor_centres=[a+angle for a in CFG['rotor_slot_angles_deg']] if level>=3 else [] circle(bore,arc_deg=arc_deg,mouths=stator_centres,half_angle=PROFILES['stator']['mouth_half_angle_deg']) circle(rotor,eccentricity,arc_deg=arc_deg,mouths=rotor_centres,half_angle=PROFILES['rotor']['mouth_half_angle_deg']) circle(shaft, eccentricity, arc_deg=arc_deg) label(outer - 8, 0, 'Воздух и немагнитный вал', 10, 4) label((stator+split)/2, 0, 'Сталь спинки статора — условная', steel_mesh, 1) label((bore+CFG['stator_slot_radii_mm'][0])/2,0,'Сталь зубцов статора — условная',steel_mesh,1) label((shaft+CFG['rotor_slot_radii_mm'][0])/2 + eccentricity, 0, 'Сталь ротора — условная', steel_mesh, 2) label(eccentricity, 0, 'Воздух и немагнитный вал', steel_mesh, 2) label(-(rotor+bore)/2, 0, 'Воздух и немагнитный вал', gap_mesh, 3) physical_currents = [-SQRT2 * current * math.cos((current_angle - k*120)*PI/180) for k in range(3)] # Ток I направлен из генератора; ток FEMM направлен в обмотку и имеет обратный знак. if id_axis is not None: physical_currents = [id_axis[0] * math.cos((angle-k*120)*PI/180) - id_axis[1] * math.sin((angle-k*120)*PI/180) for k in range(3)] for name, amps in zip(PHASES, physical_currents): mi_addcircprop(name, amps if level >= 5 else 0, 1) mi_addcircprop(FIELD, excitation if level >= 5 else 0, 1) if level >= 2: # Шесть 60-градусных зон: A+, C-, B+, A-, C+, B- по пространственному углу. # Магнитная ось фазы A — x; положительная сторона A располагается сверху. bands = [(0, 1), (2, -1), (1, 1), (0, -1), (2, 1), (1, -1)] for k in range(CFG['stator_slots']): a = (k+0.5)*360/CFG['stator_slots'] slot('stator',a) wedge_point=local_to_world([(PROFILES['stator']['xb_mm']+PROFILES['stator']['xe_mm'])/2,0],a) label(wedge_point[0],wedge_point[1],'Клин статора — немагнитный, условный',slot_mesh,13) band = int(((a - 60) % 360) / 60) phase, sign = bands[band] p = xy(sum(CFG['stator_slot_radii_mm'])/2, a) label(p[0], p[1], 'Медь обмотки' if level >= 5 else 'Воздух и немагнитный вал', slot_mesh, 10 + phase, PHASES[phase] if level >= 5 else '', sign * CFG['slot_turns_stator'] if level >= 5 else 0) if level >= 3: for base_angle in CFG['rotor_slot_angles_deg']: a = base_angle + angle slot('rotor',a,eccentricity) wedge_point=local_to_world([(PROFILES['rotor']['xb_mm']+PROFILES['rotor']['xe_mm'])/2,0],a,eccentricity) label(wedge_point[0],wedge_point[1],'Клин ротора — немагнитный, условный',slot_mesh,2) p = xy(sum(CFG['rotor_slot_radii_mm'])/2, a, eccentricity) label(p[0], p[1], 'Медь обмотки' if level >= 5 else 'Воздух и немагнитный вал', slot_mesh, 2, FIELD if level >= 5 else '', (CFG['slot_turns_field'] if base_angle < 180 else -CFG['slot_turns_field']) if level >= 5 else 0) return physical_currents def solve(tag, angle=-90, excitation=None, current=0, current_angle=0, mesh=1, eccentricity=0, nonlinear=True, level=5, outer=None, id_axis=None, scan=False): global serial if excitation is None: excitation=CFG['baseline_excitation_A'] serial += 1 amps = build(level, angle, excitation, current, current_angle, mesh, eccentricity, nonlinear, outer, id_axis) stem = '%03d-open-wedge-%s' % (serial, tag) path = OUT + '/' + stem + '.fem' mi_saveas(path) mi_close() opendocument(path) mi_analyze(1) mi_loadsolution() psi = [real(mo_getcircuitproperties(name)[2]) for name in PHASES] space = complex((2*psi[0]-psi[1]-psi[2])/3, (psi[1]-psi[2])/math.sqrt(3)) ip = complex(current * math.cos(current_angle*PI/180), current * math.sin(current_angle*PI/180)) voltage = 1j * OMEGA * space / SQRT2 - CFG['phase_resistance_ohm'] * ip power = 3 * (voltage.real * ip.real + voltage.imag * ip.imag) pf = power / (3 * abs(voltage) * current) if current > 0 and abs(voltage) > 1e-10 else 0 mo_clearblock() mo_groupselectblock(2) torque = real(mo_blockintegral(22)) fx = real(mo_blockintegral(18)) fy = real(mo_blockintegral(19)) mo_clearblock() scan_radius=(CFG['rotor_radius_mm']+CFG['bore_radius_mm'])/2 probe = mo_getpointvalues(scan_radius*math.cos(20*PI/180), scan_radius*math.sin(20*PI/180)) row = {'tag': tag, 'model': stem+'.fem', 'solution': stem+'.ans', 'geometry_revision':CFG['geometry_revision'], 'saved_reopened_solved': True, 'level': level, 'mesh_level': mesh, 'rotor_angle_deg': angle, 'excitation_A': excitation, 'current_rms_A': current, 'outgoing_current_angle_deg': current_angle, 'physical_phase_currents_A': amps, 'eccentricity_mm': eccentricity, 'nonlinear': nonlinear, 'elements': mo_numelements(), 'nodes': mo_numnodes(), 'phase_flux_Wb': psi, 'flux_space_Wb': [space.real, space.imag], 'voltage_V': [voltage.real, voltage.imag], 'voltage_rms_V': abs(voltage), 'power_factor': pf, 'terminal_power_W': power, 'copper_loss_W': 3*CFG['phase_resistance_ohm']*current*current, 'rotor_torque_Nm': torque, 'rotor_force_N': [fx, fy], 'probe_B_T': [real(probe[1]), real(probe[2])]} if tag in ['02-checkpoint','03-checkpoint','04-checkpoint','05-checkpoint']: geometry_probes=[] geometry_areas=[] for kind in (['stator','rotor'] if level>=3 else ['stator']): p=PROFILES[kind] centre=GEOMETRY_CONTRACT['stator_slot_centres_deg'][0] if kind=='stator' else CFG['rotor_slot_angles_deg'][0]+angle ecc=0 if kind=='stator' else eccentricity for region,local in p['probes_local_mm'].items(): point=local_to_world(local,centre,ecc) values=mo_getpointvalues(point[0],point[1]) mur=[real(values[9]),real(values[10])] conductivity=real(values[3]) expected_sigma=58 if region=='copper' and level>=5 else 0 passed=max(abs(v-1) for v in mur)<1e-6 and abs(conductivity-expected_sigma)<1e-6 geometry_probes.append({'kind':kind,'region':region,'point_mm':point, 'relative_permeability':mur,'conductivity_MSm':conductivity, 'expected_relative_permeability':1,'expected_conductivity_MSm':expected_sigma,'passed':passed}) for region,area_key,local in [('wedge','wedge_area_mm2',p['probes_local_mm']['wedge']), ('copper','copper_area_mm2',p['probes_local_mm']['copper'])]: point=local_to_world(local,centre,ecc) mo_clearblock() mo_selectblock(point[0],point[1]) native_mm2=real(mo_blockintegral(5))*1e6 analytic_mm2=p[area_key] error=abs(native_mm2-analytic_mm2)/analytic_mm2 geometry_areas.append({'kind':kind,'region':region,'native_area_mm2':native_mm2, 'analytic_area_mm2':analytic_mm2,'relative_error':error,'passed':error<1e-7}) mo_clearblock() row['geometry_material_probes']=geometry_probes row['geometry_area_checks']=geometry_areas check(tag+': устья, шлицы, клинья и медь без железных перемычек',all(p['passed'] for p in geometry_probes)) check(tag+': нативные площади клина и меди',all(p['passed'] for p in geometry_areas), max(p['relative_error'] for p in geometry_areas),1e-7) if scan: samples = [] for k in range(180): a = 2*PI*k/180 p = mo_getpointvalues(scan_radius*math.cos(a), scan_radius*math.sin(a)) samples.append([k*2, real(p[1])*math.cos(a)+real(p[2])*math.sin(a), -real(p[1])*math.sin(a)+real(p[2])*math.cos(a)]) row['gap_angle_Br_Bt'] = samples if tag=='05-checkpoint': # Экспорт численных треугольников: это данные МКЭ, не снимок интерфейса FEMM. mo_smooth('off') with open(OUT+'/05-field-triangles.json','w') as f: f.write('{"schema":"femm5-field-triangles-v1","unit":"mm,T","source_model":'+json.dumps(stem+'.fem')+',"triangles":[') for k in range(1,int(mo_numelements())+1): e=mo_getelement(k) points=[] for j in range(3): px,py=mo_getnode(int(e[j])) points.extend([round(real(px),6),round(real(py),6)]) p=mo_getpointvalues(real(e[3]),real(e[4])) points.append(round(math.sqrt(real(p[1])**2+real(p[2])**2),6)) if k>1: f.write(',') f.write(json.dumps(points)) f.write(']}') row['field_data']='05-field-triangles.json' RECORDS.append(row) mo_close() mi_close() print('SOLVED', stem, row['elements'], row['voltage_rms_V']) # Промежуточная квитанция сохраняет уже завершённые расчёты при остановке серии. with open(OUT + '/progress.json', 'w') as f: json.dump({'records': RECORDS, 'completed_solves': len(RECORDS)}, f) return row def csv(name, headers, rows): with open(OUT + '/' + name + '.csv', 'w') as f: f.write(','.join(headers)+'\n') for row in rows: f.write(','.join(str(v) for v in row)+'\n') def check(name, passed, value=None, limit=None): item = {'name': name, 'passed': bool(passed)} if value is not None: item['value'] = value if limit is not None: item['limit'] = limit CHECKS.append(item) def relative(a, b): return abs(a-b)/max(abs(b), 1e-12) def operating(tag, angle=-90, excitation=None, current=0, current_angle=0, mesh=1, outer=None): """Основная гармоника: усредняем Ψspace(θ)·exp(-jθ) по шести положениям. В здоровой симметричной машине фазовая перестановка продолжает 60° до 360°. Шаг 10° эквивалентен 36 отсчётам за период; гармоники 36k±1 могут алиаситься. Этап 15 проверяет основную гармонику и момент полным периодом 72 отсчёта. В дефектной машине это сокращение не применяется. """ rows=[] if excitation is None: excitation=CFG['baseline_excitation_A'] flux=0j for k in range(6): theta=k*10 r=solve(tag+'-fundamental-%02d'%k,angle+theta,excitation,current,current_angle+theta,mesh=mesh,outer=outer) rows.append(r) flux+=complex(*r['flux_space_Wb'])*complex(math.cos(theta*PI/180),-math.sin(theta*PI/180))/6 ip=complex(current*math.cos(current_angle*PI/180),current*math.sin(current_angle*PI/180)) voltage=1j*OMEGA*flux/SQRT2-CFG['phase_resistance_ohm']*ip power=3*(voltage.real*ip.real+voltage.imag*ip.imag) result=dict(rows[0]) result['fundamental_sample_models']=[r['model'] for r in rows] result['fundamental_sample_count_full_period_equivalent']=36 result['flux_space_Wb']=[flux.real,flux.imag] result['voltage_V']=[voltage.real,voltage.imag] result['voltage_rms_V']=abs(voltage) result['terminal_power_W']=power result['power_factor']=power/(3*abs(voltage)*current) if current and abs(voltage)>1e-10 else 0 result['mean_torque_Nm']=sum(r['rotor_torque_Nm'] for r in rows)/6 result['voltage_definition']='Основная гармоника шести полевых положений за60°; эквивалент36 за360° для здоровой симметричной машины.' return result def spectrum(samples, maximum=10): """Односторонняя ДПФ: a_h=2/N Σx_k cos(hθk), b_h=2/N Σx_k sin(hθk). Амплитуда sqrt(a_h²+b_h²); E_h,rms=hωΨ_h/√2. Постоянная и Найквист исключены. """ n = len(samples) rows = [] for h in range(1, maximum+1): a = 2*sum(samples[k]*math.cos(2*PI*h*k/n) for k in range(n))/n b = 2*sum(samples[k]*math.sin(2*PI*h*k/n) for k in range(n))/n rows.append({'harmonic': h, 'cos': a, 'sin': b, 'amplitude': math.sqrt(a*a+b*b)}) return rows def target(tag, amps, pf=None, voltage=None, guess=None): """Полевой подбор If и угла к основной гармонике U=Uзад+j0. F=(Re(U)−Uзад,Im(U)); численный якобиан по приращениям If и α. Решаем J·Δx=−F как систему2×2 по определителю, ограничивая длину шага. Невязка ||F||/Uзад измеряет одновременно модуль и фазу напряжения. """ if pf is None: pf=CFG['target_power_factor'] if voltage is None: voltage=CFG['target_phase_voltage_rms_V'] beta = -math.acos(pf)*180/PI initial_current=complex(amps*math.cos(beta*PI/180),amps*math.sin(beta*PI/180)) initial_emf=complex(voltage,0)+complex(CFG['phase_resistance_ohm'],CFG['iteration_initial_reactance_ohm'])*initial_current seed=(max(.1,abs(initial_emf)*CFG['iteration_initial_excitation_A_per_V']),math.atan2(initial_emf.imag,initial_emf.real)*180/PI-90) # Грубая реактивность служит только начальным приближением; все ответы подтверждает поле. x, a = guess if guess else seed trace = [] result = None jacobian=None previous_residual=1e100 for iteration in range(12): result = operating(tag+'-it%d' % iteration, a, x, amps, beta) v = result['voltage_V'] residual = math.sqrt((v[0]-voltage)**2+v[1]**2)/voltage trace.append([iteration, x, a, v[0], v[1], residual]) if residual < 0.005: break if jacobian is None or iteration%3==0 or residual>previous_residual*.7: dx, da = max(0.5, x*0.03), 2.0 r1 = operating(tag+'-If-difference', a, x+dx, amps, beta) r2 = operating(tag+'-angle-difference', a+da, x, amps, beta) jacobian=[(r1['voltage_V'][0]-v[0])/dx,(r2['voltage_V'][0]-v[0])/da, (r1['voltage_V'][1]-v[1])/dx,(r2['voltage_V'][1]-v[1])/da] j11,j12,j21,j22=jacobian previous_residual=residual det = j11*j22-j12*j21 if abs(det) < 1e-8: raise RuntimeError('Вырожденный якобиан полевого подбора') ex, ey = voltage-v[0], -v[1] step_x = (ex*j22-j12*ey)/det step_a = (j11*ey-ex*j21)/det x = max(0.1, min(150, x+max(-15, min(15, step_x)))) a += max(-25, min(25, step_a)) csv(tag+'-iterations', ['iteration','excitation_A','angle_deg','Ureal_V','Uimag_V','residual'], trace) check(tag+': напряжение и фаза', residual < 0.005, residual, 0.005) if amps: check(tag+': коэффициент мощности', abs(result['power_factor']-pf) < 0.01, abs(result['power_factor']-pf), 0.01) result['target_relative_residual'] = residual return result def waveform(tag, excitation, current=0, base_angle=-90, beta=0, eccentricity=0, count=24): rows = [] for k in range(count): theta = k*360/count row = solve(tag+'-p%02d' % k, base_angle+theta, excitation, current, beta+theta, eccentricity=eccentricity) rows.append(row) # Конечные разности по замкнутому периоду; E — производная Ψ с выбранным знаком. dt = 1/(FREQ*count) data = [] for k, row in enumerate(rows): emf = -(rows[(k+1)%count]['phase_flux_Wb'][0]-rows[(k-1)%count]['phase_flux_Wb'][0])/(2*dt) data.append([k*360/count, k*dt, row['phase_flux_Wb'][0], emf, row['rotor_torque_Nm'], row['probe_B_T'][0], row['probe_B_T'][1], row['rotor_force_N'][0], row['rotor_force_N'][1]]) csv(tag, ['angle_deg','time_s','Psi_A_Wb','e_A_central_difference_V','torque_Nm', 'probe_Bx_T','probe_By_T','force_x_N','force_y_N'], data) harmonics = spectrum([r['phase_flux_Wb'][0] for r in rows], count//2-1) for h in harmonics: h['emf_rms_V'] = OMEGA*h['harmonic']*h['amplitude']/SQRT2 csv(tag+'-harmonics', ['harmonic','flux_amplitude_Wb','emf_rms_V'], [[h['harmonic'],h['amplitude'],h['emf_rms_V']] for h in harmonics]) return {'samples': rows, 'harmonics': harmonics, 'sample_count': count, 'period_mechanical_deg': 360, 'period_s': 1/FREQ, 'sampling_limit': 'Гармоники выше Найквиста не определены; спектр может содержать алиасинг пазовых гармоник.'} def dq_parameters(): # Park при α=0: Ψd=Re(Ψspace), Ψq=Im(Ψspace), токи — амплитудные. # Lkj=∂Ψk/∂ij≈[Ψk(+Δij)−Ψk(−Δij)]/(2Δij) при фиксированных If,α. # Секущая [Ψk(10А)−Ψk(0)]/10А отличается от производной в насыщенной стали. bias=CFG['baseline_excitation_A'] base = solve('09-dq-bias', 0, bias, id_axis=[0, 0]) def transformed(row): return row['flux_space_Wb'] origin = transformed(base) matrix = [] secants = [] for axis in range(2): plus = [0, 0] minus = [0, 0] plus[axis], minus[axis] = 1, -1 p = solve('09-dq-plus%d' % axis, 0, bias, id_axis=plus) m = solve('09-dq-minus%d' % axis, 0, bias, id_axis=minus) matrix.append([(transformed(p)[k]-transformed(m)[k])/2 for k in range(2)]) big = [0, 0] big[axis] = 10 r = solve('09-dq-secant%d' % axis, 0, bias, id_axis=big) secants.append((transformed(r)[axis]-origin[axis])/10) # matrix[возмущаемый ток][измеряемое потокосцепление], единицы Гн. reciprocal = abs(matrix[0][1]-matrix[1][0])/max(abs(matrix[0][0]),abs(matrix[1][1]),1e-12) check('09: взаимность дифференциальной матрицы', reciprocal < 0.05, reciprocal, 0.05) check('09: положительные Ldd и Lqq', matrix[0][0]>0 and matrix[1][1]>0) data = {'bias_field_A':bias, 'increment_peak_A':1, 'secant_peak_A':10, 'field_flux_d_Wb':origin[0], 'field_flux_q_Wb':origin[1], 'differential_H': {'Ldd':matrix[0][0], 'Ldq':matrix[1][0], 'Lqd':matrix[0][1], 'Lqq':matrix[1][1]}, 'secant_H': {'Ld':secants[0], 'Lq':secants[1]}, 'Xd_differential_ohm': OMEGA*matrix[0][0], 'Xq_differential_ohm': OMEGA*matrix[1][1], 'reciprocity_relative_error':reciprocal, 'scope':'Локальные параметры 2D активной части при If=30 А, фиксированный ротор; без лобового рассеяния.'} csv('09-dq', ['Ldd_H','Ldq_H','Lqd_H','Lqq_H','Ld_secant_H','Lq_secant_H'], [[matrix[0][0],matrix[1][0],matrix[0][1],matrix[1][1],secants[0],secants[1]]]) return data def reduced_dq(parameters): """Связанные электрические и механические уравнения автономного генератора. Ψd=Ψf+Ld·id, Ψq=Lq·iq; vd=R·id+Ld·did/dt−Ω·Ψq, vq=R·iq+Lq·diq/dt+Ω·Ψd. Пассивные токи id,iq направлены в машину. Для симметричной резистивной нагрузки v=−Rн·i. Момент двигателя Te=3p/2·(Ψd·iq−Ψq·id), поэтому у генератора Te<0. J·dΩ/dt=Tm+Te−DΩ; dδ/dt=Ω−ωs. Число пар полюсов p=1. Это постоянные локальные параметры активной части, не новая динамическая МКЭ. """ ld=parameters['differential_H']['Ldd'] lq=parameters['differential_H']['Lqq'] flux=parameters['field_flux_d_Wb'] resistance=CFG['phase_resistance_ohm'] inertia=.2 damping=.001 initial_load=8.0 total=resistance+initial_load iq0=-OMEGA*flux/(total+OMEGA*OMEGA*ld*lq/total) id0=OMEGA*lq*iq0/total def torque(state): d,q=state[0],state[1] return 1.5*((flux+ld*d)*q-lq*q*d) initial=[id0,iq0,OMEGA,0.0] tm=-torque(initial)+damping*OMEGA def rhs(t,state): d,q,speed,delta=state load=8.0 if t<.1 else 6.0 rd=resistance+load return [(-rd*d+speed*lq*q)/ld, (-rd*q-speed*(ld*d+flux))/lq, (tm+torque(state)-damping*speed)/inertia, speed-OMEGA] def integrate(dt): state=list(initial) result=[] n=int(.5/dt+.5) for k in range(n+1): t=k*dt load=8.0 if t<.1 else 6.0 power=1.5*load*(state[0]**2+state[1]**2) result.append([t,state[0],state[1],state[2],state[3]*180/PI, torque(state),tm,power,load]) if k==n: break # RK4: четыре согласованных вычисления правой части за шаг. a=rhs(t,state) b=rhs(t+dt/2,[state[j]+dt*a[j]/2 for j in range(4)]) c=rhs(t+dt/2,[state[j]+dt*b[j]/2 for j in range(4)]) # В точке скачка используем левый предел для конца предыдущего шага. end=t+dt-1e-12 if t<.1 and t+dt>=.1 else t+dt d=rhs(end,[state[j]+dt*c[j] for j in range(4)]) state=[state[j]+dt*(a[j]+2*b[j]+2*c[j]+d[j])/6 for j in range(4)] return result coarse=integrate(.0002) fine=integrate(.0001) mismatch=max(abs(coarse[k][j]-fine[k*2][j])/max(abs(fine[k*2][j]),1.0) for k in range(len(coarse)) for j in [1,2,3]) initial_rhs=rhs(0,initial) steady=max(abs(v) for v in initial_rhs) check('16: начальное электрическое и механическое равновесие',steady<1e-8,steady,1e-8) check('16: уточнение RK4 dt/2, токи и скорость',mismatch<.005,mismatch,.005) check('16: генераторный момент и положительная мощность нагрузки', all(r[5]<0 and r[7]>0 and r[3]>0 for r in fine)) headers=['time_s','id_peak_A','iq_peak_A','speed_rad_s','delta_deg','Te_Nm','Tm_Nm','load_power_W','load_ohm'] csv('16-dq-dynamics',headers,fine) csv('16-dq-dynamics-coarse',headers,coarse) static_comparison=[] for k in [0,len(fine)//2,len(fine)-1]: sample=fine[k] native=solve('16-dq-validation-%04d'%k,angle=0,excitation=CFG['baseline_excitation_A'],id_axis=[sample[1],sample[2]]) predicted=[flux+ld*sample[1],lq*sample[2]] actual=native['flux_space_Wb'] error=abs(complex(*predicted)-complex(*actual))/max(abs(complex(*actual)),1e-12) static_comparison.append({'time_s':sample[0],'id_peak_A':sample[1],'iq_peak_A':sample[2], 'frozen_flux_dq_Wb':predicted,'native_flux_dq_Wb':actual, 'relative_flux_discrepancy':error,'record':native}) csv('16-dq-static-comparison',['time_s','id_peak_A','iq_peak_A','frozen_Psid_Wb','frozen_Psiq_Wb', 'native_Psid_Wb','native_Psiq_Wb','relative_flux_discrepancy'], [[r['time_s'],r['id_peak_A'],r['iq_peak_A']]+r['frozen_flux_dq_Wb']+r['native_flux_dq_Wb']+[r['relative_flux_discrepancy']] for r in static_comparison]) return {'Ld_H':ld,'Lq_H':lq,'field_flux_d_Wb':flux,'R_ohm':resistance, 'J_kg_m2':inertia,'D_Nm_s':damping,'mechanical_torque_Nm':tm, 'initial_state_id_iq_speed_delta_rad':initial,'initial_rhs_max':steady, 'load_step_ohm':[8,6],'load_step_time_s':.1,'time_step_s':.0001, 'dt_refinement_relative_error':mismatch,'final_state':fine[-1], 'static_nonlinear_comparison':static_comparison, 'trajectory_current_range_peak_A':{'id':[min(r[1] for r in fine),max(r[1] for r in fine)], 'iq':[min(r[2] for r in fine),max(r[2] for r in fine)]}, 'nonlinear_transient_validated':False, 'static_comparison_scope':'Три нелинейных магнитостатических сравнения Ψ(i) в фиксированном угловом положении; не проверка всей переходной траектории и не физический эксперимент.', 'parameter_calibration_warning':'Ld,Lq линеаризованы центральной разностью ±1 А около id=iq=0; фактические диапазоны траектории приведены в trajectory_current_range_peak_A. Это экстраполяция локальных параметров. Три статических нелинейных МКЭ сравнения отдельно количественно показывают ошибку замороженного отображения Ψ(i); траектория не является подтверждённым нелинейным переходным процессом.', 'secant_vs_differential_relative_difference':{'d':relative(parameters['secant_H']['Ld'],ld), 'q':relative(parameters['secant_H']['Lq'],lq)}, 'scope':'Автономный генератор на симметричную резистивную нагрузку. Замороженные локальные Ld,Lq,Ψf из статических МКЭ; без демпфера, динамики возбуждения, АРВ, сети и переменного насыщения. Инерция и демпфирование условные.'} def stage(stage_id): if stage_id in ['01','02','03','04','05']: level = int(stage_id) row = solve(stage_id+'-checkpoint', level=level, scan=level==5) if level==5: check('05: ненулевое потокосцепление возбуждённой машины', abs(complex(*row['flux_space_Wb']))>1e-5) return row if stage_id == '06': rows = [operating('06-mesh%d' % k, mesh=k) for k in range(3)] delta = relative(rows[1]['voltage_rms_V'],rows[2]['voltage_rms_V']) check('06: последние две сетки, напряжение ХХ', delta<0.02,delta,0.02) outer = operating('06-outer-boundary',mesh=2,outer=CFG['boundary_radius_mm']+30) drift = relative(outer['voltage_rms_V'],rows[2]['voltage_rms_V']) check('06: удаление внешней границы',drift<0.01,drift,0.01) csv('06-mesh',['mesh_level','elements','voltage_rms_V'], [[r['mesh_level'],r['elements'],r['voltage_rms_V']] for r in rows]) return {'meshes':rows,'outer_boundary':outer,'last_mesh_relative_change':delta,'boundary_relative_change':drift, 'scope':'Сходимость только интегрального напряжения ХХ данного режима; не всех сил/моментов и насыщенных режимов.'} if stage_id == '07': wave = waveform('07-no-load',CFG['baseline_excitation_A'],count=36 if CASE=='standard' else 12) sample = solve('07-gap-scan',scan=True) csv('07-gap',['angle_deg','Br_T','Bt_T'],sample['gap_angle_Br_Bt']) wave['gap_spectrum'] = spectrum([p[1] for p in sample['gap_angle_Br_Bt']],15) check('07: основная гармоника потокосцепления',wave['harmonics'][0]['amplitude']>0.01) return wave if stage_id == '08': return target('08-load',CFG['target_current_rms_A']) if stage_id == '09': return dq_parameters() if stage_id == '10': rows = [operating('10-occ-%02d' % n,excitation=n) for n in CFG['occ_excitation_A']] csv('10-occ',['excitation_A','voltage_rms_V','flux_A_Wb'], [[r['excitation_A'],r['voltage_rms_V'],r['phase_flux_Wb'][0]] for r in rows]) check('10: монотонность ХХ', all(a['voltage_rms_V']1e-12 else x1 last = operating('11-scc-If%02d-it%d' % (exc,it),excitation=exc,current=x2,current_angle=-90) y2 = last['flux_space_Wb'][1] if abs(y2) < 0.0001: break x0,y0,x1,y1=x1,y1,x2,y2 residual = abs(complex(*last['flux_space_Wb']))/max(abs(complex(*a['flux_space_Wb'])),1e-12) check('11: КЗ If=%d А, остаточный поток' % exc,residual<0.01,residual,0.01) rows.append({'excitation_A':exc,'current_rms_A':last['current_rms_A'], 'flux_residual_relative':residual,'record':last}) csv('11-scc',['excitation_A','current_rms_A','flux_residual_relative'], [[r['excitation_A'],r['current_rms_A'],r['flux_residual_relative']] for r in rows]) return {'samples':rows,'resistance_ohm_for_scc_constraint':0, 'scope':'Симметричное установившееся КЗ с R=0 в условии Ψd=Ψq≈0. Не ударный/переходный ток. Напряжение с R=0.15 из общей записи не служит критерием КЗ.'} if stage_id == '12': rows=[] # При постоянных If и относительном угле тока к ЭДС получаем нагрузочную кривую. for amps in CFG['load_currents_A']: rows.append(operating('12-load-%02d' % amps,current=amps,current_angle=-math.acos(CFG['target_power_factor'])*180/PI)) csv('12-load',['current_rms_A','voltage_rms_V','actual_power_factor','terminal_power_W'], [[r['current_rms_A'],r['voltage_rms_V'],r['power_factor'],r['terminal_power_W']] for r in rows]) return {'samples':rows,'scope':'Постоянные If и угол тока относительно ЭДС ХХ; фактический cosφ к напряжению изменяется и записан.'} if stage_id == '13': rows=[] guess=None for amps in CFG['regulation_currents_A']: r=RESULTS['08'] if amps==CFG['target_current_rms_A'] and '08' in RESULTS else target('13-regulation-%02d' % amps,amps,guess=guess) guess=(r['excitation_A'],r['rotor_angle_deg']) rows.append(r) csv('13-regulation',['current_rms_A','excitation_A','rotor_angle_deg','voltage_rms_V','power_factor','target_residual'], [[r['current_rms_A'],r['excitation_A'],r['rotor_angle_deg'],r['voltage_rms_V'],r['power_factor'],r['target_relative_residual']] for r in rows]) return {'samples':rows,'target_voltage_rms_V':CFG['target_phase_voltage_rms_V'],'target_power_factor':CFG['target_power_factor'],'field_constrained':True} if stage_id == '14': # Для заданных U и P перебираем Q: I=(P-jQ)/(3U), каждый узел подтверждаем полем. rows=[] base_power=3*CFG['target_phase_voltage_rms_V']*CFG['target_current_rms_A'] for power in [base_power*.5,base_power]: guess=None for reactive in [-base_power*.5,0,base_power*.5,base_power]: active_i=power/(3*CFG['target_phase_voltage_rms_V']) reactive_i=reactive/(3*CFG['target_phase_voltage_rms_V']) current=math.sqrt(active_i**2+reactive_i**2) pf=active_i/current beta=-math.atan2(reactive_i,active_i)*180/PI # target принимает отстающий cosφ; опережающий ток задаём отдельным знаком. if reactive < 0: r=target_signed('14-U-P%d-Q%d' % (power,reactive),current,beta,guess) else: r=target('14-U-P%d-Q%d' % (power,reactive),current,pf,guess=guess) guess=(r['excitation_A'],r['rotor_angle_deg']) rows.append({'target_P_W':power,'target_Q_var':reactive,'record':r}) csv('14-u-curves',['target_P_W','target_Q_var','excitation_A','current_rms_A','actual_P_W','voltage_rms_V','residual'], [[r['target_P_W'],r['target_Q_var'],r['record']['excitation_A'],r['record']['current_rms_A'], r['record']['terminal_power_W'],r['record']['voltage_rms_V'],r['record']['target_relative_residual']] for r in rows]) return {'samples':rows,'target_voltage_rms_V':CFG['target_phase_voltage_rms_V'],'field_constrained':True, 'scope':'Редкое учебное семейство по 4 точки; в узлах U и P заданы через полевой подбор. Для точного минимума необходима более частая сетка Q.'} if stage_id == '15': r=RESULTS['08'] if '08' in RESULTS else target('15-base-load',CFG['target_current_rms_A']) current=CFG['target_current_rms_A'] beta=-math.acos(CFG['target_power_factor'])*180/PI count=72 if CASE=='standard' else 36 a=waveform('15-healthy',r['excitation_A'],current,r['rotor_angle_deg'],beta,count=count) b=waveform('15-eccentric',r['excitation_A'],current,r['rotor_angle_deg'],beta,CFG['eccentricity_mm'],count) # В полном периоде сохраняется статическое смещение оси; центр ротора не вращается. for name,wave in [('healthy',a),('eccentric',b)]: torque=sum(s['rotor_torque_Nm'] for s in wave['samples'])/count flux_h=wave['harmonics'] wave['mean_torque_Nm']=torque wave['even_flux_ratio']=math.sqrt(sum(h['amplitude']**2 for h in flux_h if h['harmonic']%2==0))/max(flux_h[0]['amplitude'],1e-12) check('15: минимальный зазор положителен',CFG['bore_radius_mm']-CFG['rotor_radius_mm']-CFG['eccentricity_mm']>0) # Проверка энергии использует средний электромагнитный момент и активную мощность. expected=r['terminal_power_W']+r['copper_loss_W'] mechanical=-a['mean_torque_Nm']*OMEGA balance=relative(mechanical,expected) check('15: средний момент и электрическая мощность',balance<0.12,balance,0.12) harmonic=a['harmonics'][0] lambda1=complex(harmonic['cos'],-harmonic['sin']) iout=complex(current*math.cos(beta*PI/180),current*math.sin(beta*PI/180)) full_voltage=1j*OMEGA*lambda1/SQRT2-CFG['phase_resistance_ohm']*iout full_power=3*(full_voltage.real*iout.real+full_voltage.imag*iout.imag) full_residual=abs(full_voltage-complex(CFG['target_phase_voltage_rms_V'],0))/CFG['target_phase_voltage_rms_V'] coarse=spectrum([s['phase_flux_Wb'][0] for s in a['samples'][::2]],10)[0] voltage_refinement=relative(harmonic['amplitude'],coarse['amplitude']) torque_coarse=sum(s['rotor_torque_Nm'] for s in a['samples'][::2])/(count/2) torque_refinement=relative(a['mean_torque_Nm'],torque_coarse) check('15: основная гармоника полного периода к заданному U',full_residual<.03,full_residual,.03) check('15: уточнение угловой сетки, основная гармоника',voltage_refinement<.02,voltage_refinement,.02) check('15: уточнение угловой сетки, средний момент',torque_refinement<.03,torque_refinement,.03) return {'healthy':a,'eccentric':b,'base_load':r,'mechanical_power_from_mean_torque_W':mechanical, 'electromagnetic_power_from_flux_W':expected,'power_balance_relative_error':balance, 'full_period_voltage_V':[full_voltage.real,full_voltage.imag], 'full_period_power_W':full_power,'full_period_target_residual':full_residual, 'full_period_power_balance_relative_error':relative(mechanical,full_power+r['copper_loss_W']), 'sampling_refinement':{'coarse_samples':count//2,'fine_samples':count, 'fundamental_relative_change':voltage_refinement,'mean_torque_relative_change':torque_refinement}, 'scope':'Квазистатическая последовательность; статический эксцентриситет 0.4 мм. Диагностический признак относится только к этому опыту, не универсальный классификатор дефекта.'} if stage_id == '16': # Форма P(δ) получена из поля при фиксированных токах; инерция и демпфирование условные. rows=[operating('16-power-angle-%03d' % a,angle=-180+a,current=CFG['target_current_rms_A'],current_angle=0) for a in [0,15,30,45,60,75,90]] table=[[r['rotor_angle_deg']+180,r['terminal_power_W']+r['copper_loss_W']] for r in rows] # Только редуцированная демонстрация: полевая таблица, кусочно-линейная интерполяция. def pfield(delta): d=delta*180/PI for j in range(len(table)-1): if table[j][0] <= d <= table[j+1][0]: t=(d-table[j][0])/(table[j+1][0]-table[j][0]) return table[j][1]*(1-t)+table[j+1][1]*t return table[0][1] if d=0.1 else 1) pe=pfield(delta) series.append([t,delta*180/PI,speed,input_power,pe]) if k<1000: acceleration=(input_power-pe)/(inertia*OMEGA)-damping*speed/inertia speed+=dt*acceleration delta+=dt*speed csv('16-power-angle',['angle_deg','electromagnetic_power_W'],table) csv('16-swing',['time_s','delta_deg','delta_dot_rad_s','Pm_W','Pe_W'],series) check('16: траектория в диапазоне полевой таблицы',all(0<=s[1]<=90 for s in series)) dq=RESULTS.get('09',None) if dq is None: dq=dq_parameters() coupled=reduced_dq(dq) return {'coupled_dq_model':coupled,'field_power_angle':table,'inertia_kg_m2':inertia,'damping_Nm_s':damping, 'time_step_s':dt,'final_delta_deg':series[-1][1], 'scope':'Редуцированное уравнение J·ωs·δ̈+D·ωs·δ̇=Pm−Pe(δ), фиксированные статорные токи. Не динамика подключённого к сети ТГ: нет электрических дифференциальных уравнений, АРВ, демпфера, насыщения во времени и замыкания сети.'} def target_signed(tag, amps, beta, guess=None): """То же полевое условие U=80+j0 для опережающего тока.""" # Смена временного знака угла тока не меняет построение модели или определения V. voltage=CFG['target_phase_voltage_rms_V'] initial_current=complex(amps*math.cos(beta*PI/180),amps*math.sin(beta*PI/180)) initial_emf=complex(voltage,0)+complex(CFG['phase_resistance_ohm'],CFG['iteration_initial_reactance_ohm'])*initial_current seed=(max(.1,abs(initial_emf)*CFG['iteration_initial_excitation_A_per_V']),math.atan2(initial_emf.imag,initial_emf.real)*180/PI-90) x,a=guess if guess else seed trace=[] for it in range(12): r=operating(tag+'-it%d'%it,a,x,amps,beta) vx,vy=r['voltage_V'] residual=math.sqrt((vx-voltage)**2+vy**2)/voltage trace.append([it,x,a,vx,vy,residual]) if residual<0.005: break dx=max(.5,x*.03) p=operating(tag+'-If-difference',a,x+dx,amps,beta) q=operating(tag+'-angle-difference',a+2,x,amps,beta) aa=(p['voltage_V'][0]-vx)/dx cc=(p['voltage_V'][1]-vy)/dx bb=(q['voltage_V'][0]-vx)/2 dd=(q['voltage_V'][1]-vy)/2 determinant=aa*dd-bb*cc if abs(determinant)<1e-8: raise RuntimeError('Вырожденный якобиан полевого подбора опережающего тока') x=max(.1,min(150,x+max(-15,min(15,((voltage-vx)*dd+bb*vy)/determinant)))) a+=max(-25,min(25,(-aa*vy-cc*(voltage-vx))/determinant)) csv(tag+'-iterations',['iteration','excitation_A','angle_deg','Ureal_V','Uimag_V','residual'],trace) r['target_relative_residual']=residual check(tag+': напряжение и фаза',residual<.005,residual,.005) return r selected = [s[0] for s in STAGES] if STAGE=='all' else ['%02d' % int(STAGE)] if any(s not in [p[0] for p in STAGES] for s in selected): raise ValueError('STAGE: all или 01..16') for stage_id in selected: RESULTS[stage_id]=stage(stage_id) with open(OUT+'/stage-'+stage_id+'.json','w') as f: json.dump(RESULTS[stage_id],f) check('Все полевые отсчёты конечны',all(math.isfinite(r['voltage_rms_V']) and math.isfinite(r['rotor_torque_Nm']) for r in RECORDS)) result={'schema':'femm5-turbogenerator-teaching-v2','geometry_revision':CFG['geometry_revision'],'runtime':sys.version, 'implementation':str(sys.implementation),'stages':STAGES,'selected_stages':selected, 'case':CASE,'parameters':CFG,'geometry_contract':GEOMETRY_CONTRACT,'records':RECORDS,'results':RESULTS,'checks':CHECKS, 'passed':all(c['passed'] for c in CHECKS),'completed_solves':len(RECORDS), 'origin':'Самостоятельные учебные размеры и результаты. Профиль открытого статора адаптирован по опубликованному Lua Милыха и Поляковой (2015, с.11); ротор зеркально адаптирован, отдельная граница воздух/клин авторская. Полный исходный Lua и исходные размеры реального ТГ не предоставлены.', 'scope':'Нелинейная 2D магнитостатика активной части и отдельная редуцированная динамика. Нет физической/3D/промышленной валидации. Числа условного образца не являются параметрами реального ТГ.'} with open(OUT+'/result.json','w') as f: json.dump(result,f) print('RESULT',OUT,'PASS' if result['passed'] else 'FAIL',len(RECORDS)) if not result['passed']: raise RuntimeError('Проверки не пройдены; изучите result.json. Не изменяйте допуски ради PASS.')