Статический эксцентриситет

"""Полные циклы несимметрии и статического эксцентриситета, без тепловой модели."""
import sys,json,math
ROOT=globals().get('TVV_ROOT','')
if ROOT in sys.path:sys.path.remove(ROOT)
sys.path.insert(0,ROOT)
for module in ('build_study','build_series','study_physics'):sys.modules.pop(module,None)
from build_study import build
from study_physics import sequence_currents,torque_about_rotor,spatial_components,cycle_summary
from physics import abc_to_dq0
from run_record import start,finish,save_json,sha
family=globals().get('TVV_STUDY','balanced')
names=['run_studies.py','build_study.py','build_series.py','geometry.py','physics.py','study_physics.py','run_record.py','machine.json','winding.json','rotor-bh.json','stator-bh.json','study-plan.json','model-scope.json']
resume=globals().get('TVV_RESUME_DIR','')
if resume:
    folder=resume
    with open(folder+'/receipt.json') as f:record=json.loads(f.read())
    if record['family']!=family:raise ValueError('Другая серия')
    for name,digest in record['inputs_and_programs'].items():
        if sha(ROOT+'/'+name)!=digest or sha(folder+'/snapshot/'+name)!=digest:raise ValueError('Изменены входы')
else:
    folder,record=start(ROOT,globals().get('TVV_OUTPUT_ROOT',ROOT+'/runs'),'study',names,'full angular cycles; imposed currents; no rotor heating claim')
    record.update(family=family,points={})
def load(name):
    with open(folder+'/snapshot/'+name+'.json') as f:return json.loads(f.read())
plan=load('study-plan');case=plan['cases'][family];machine=load('machine');winding=load('winding')['rows']
materials={'rotor':load('rotor-bh'),'stator':load('stator-bh')};rows=[]
try:
    for k in range(case['angles']):
        stem='angle-%03d'%k
        if stem in record['points']:
            for name,digest in record['points'][stem].items():
                if sha(folder+'/'+name)!=digest:raise ValueError('Повреждена точка серии')
            with open(folder+'/'+stem+'.json') as f:rows.append(json.loads(f.read()))
            continue
        theta=2*math.pi*k/case['angles'];degrees=360*k/case['angles']
        phase=sequence_currents(theta,case['positive_peak_A'],case['negative_peak_A'],case['positive_phase_rad'],case['negative_phase_rad'])
        if abs(sum(phase))>1e-7:raise ValueError('Токи не согласованы с изолированной нейтралью')
        dq=abc_to_dq0(phase,theta);ecc=case['eccentricity_mm'];path=folder+'/'+stem+'.fem'
        build(globals(),machine,winding,materials,path,case['mesh'],degrees,[dq[0],dq[1],case['If_A']],ecc)
        mi_close();opendocument(path);mi_analyze(1);mi_loadsolution();mo_smooth('off')
        names_c=[p+str(b)+' — ветвь статора' for p in 'ABC' for b in (1,2)]+['F — возбуждение ротора']
        circuits=[list(mo_getcircuitproperties(n)) for n in names_c]
        expected=[phase[j//2]/2 for j in range(6)]+[case['If_A']]
        if max(abs(a[0]-b) for a,b in zip(circuits,expected))>1e-7:raise ValueError('Неверный ток цепи')
        fluxabc=[(circuits[2*j][2]+circuits[2*j+1][2])/2 for j in range(3)]
        mo_groupselectblock();energy=mo_blockintegral(2);coenergy=mo_blockintegral(17);mo_clearblock()
        identity=abs(energy+coenergy-sum(v[0]*v[2] for v in circuits))/max(1.,abs(energy+coenergy))
        if identity>1e-5:raise ValueError('Энергетическое тождество')
        mo_groupselectblock(10);mo_groupselectblock(11);fx=mo_blockintegral(18);fy=mo_blockintegral(19);torque=mo_blockintegral(22);mo_clearblock()
        # Same air circle, two explicitly different observation frames.
        radius=machine['geometry']['rotor_body_diameter_mm']/2+machine['geometry']['air_gap_mm']/2
        samples=[]
        for j in range(plan['spatial_samples']):
            alpha=2*math.pi*j/plan['spatial_samples'];field=mo_getpointvalues(radius*math.cos(alpha),radius*math.sin(alpha))
            shifted=alpha+theta;rotorfield=mo_getpointvalues(ecc+radius*math.cos(shifted),radius*math.sin(shifted))
            samples.append({'alpha_deg':360*j/plan['spatial_samples'],'stator_BrBt_T':spatial_components(field[1],field[2],alpha),
                            'rotating_observer_BrBt_T':spatial_components(rotorfield[1],rotorfield[2],shifted)})
        row={'theta_deg':degrees,'currents_abc_A':phase,'currents_dqf_A':[dq[0],dq[1],case['If_A']],
             'lambda_abc_Wb':fluxabc,'lambda_dq0_Wb':abc_to_dq0(fluxabc,theta),'lambda_f_Wb':circuits[-1][2],
             'Fx_N':fx,'Fy_N':fy,'torque_origin_Nm':torque,'torque_rotor_Nm':torque_about_rotor(torque,fx,fy,ecc/1000),
             'energy_J':energy,'coenergy_J':coenergy,'energy_identity_error':identity,'samples':samples,'elements':mo_numelements()}
        mo_close();mi_close();save_json(folder+'/'+stem+'.json',row);rows.append(row)
        record['points'][stem]={stem+ext:sha(folder+'/'+stem+ext) for ext in ('.json','.fem','.ans')}
        save_json(folder+'/receipt.json',record);save_json(folder+'/results.json',rows);print('TVV_STUDY_POINT',family,k)
    full=cycle_summary(rows);halfrows=[dict(r,theta_deg=360*j/(len(rows)//2)) for j,r in enumerate(rows[::2])];half=cycle_summary(halfrows)
    checks={key:{stat:abs(full[key][stat]-half[key][stat])/max(abs(full[key][stat]),plan['absolute_scale']) for stat in full[key]} for key in full}
    save_json(folder+'/angular-check.json',{'full':full,'half':half,'scaled_errors':checks})
    if max(v for item in checks.values() for v in item.values())>plan['force_torque_relative']:
        raise ValueError('Средние/пульсации не сошлись по углу: нужна более густая полная серия')
    finish(folder,record,'passed_full_cycle_and_angular_check_pending_mesh_comparison')
except Exception as error:record['error']=str(error);finish(folder,record,'failed');raise
print('TVV_STUDY_RESULT',folder)