"""Учебное исследование магнитного поля во встроенном Python FEMM 5. Геометрия — мм, результаты — СИ, гармонические токи — амплитуды. Каждый запуск создаёт calculated-NN; прежние модели и решения сохраняются. Каждая модель сохраняется, закрывается, повторно открывается и решается. Группы: 0 — воздух, 1 — сердечник, 2 — исследуемое тело, 3/4 — обмотки. Условная B(H) взята из учебного комплекта курса, это не марка стали. Источники методики и точные допущения приведены в Порядок-работы.md. """ import os import sys import math import cmath import json MU0 = 4e-7 * math.pi BH = [(0, 0), (.3, 100), (.6, 200), (1, 350), (1.3, 700), (1.5, 1800), (1.65, 5000), (1.8, 15000), (2, 60000), (2.2, 200000)] source = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '') HERE = os.path.dirname(os.path.abspath(source)) if source else os.getcwd() counter = 1 while os.path.exists(HERE + '/calculated-%02d' % counter): counter += 1 OUT = HERE + '/calculated-%02d' % counter os.mkdir(OUT) ROWS = [] CHECKS = [] NODES = {} EDGES = {} def real(value): return complex(value).real def parts(value): return [complex(value).real, complex(value).imag] def rel(a, b, floor=1e-12): return abs(a - b) / max(abs(a), abs(b), floor) def check(name, error, limit): CHECKS.append({'name': name, 'error': error, 'limit': limit, 'passed': math.isfinite(error) and error <= limit}) def edge(a, b, size=0, boundary=''): for point in (a, b): if point not in NODES: mi_addnode(*point) NODES[point] = True key = tuple(sorted((a, b))) if key not in EDGES: mi_addsegment(a[0], a[1], b[0], b[1]) EDGES[key] = True mi_selectsegment((a[0] + b[0]) / 2, (a[1] + b[1]) / 2) mi_setsegmentprop(boundary, size, 0 if size else 1, 0, 0) mi_clearselected() def polygon(points, size=0): for k in range(len(points)): edge(points[k], points[(k + 1) % len(points)], size) def rect(x0, y0, x1, y1, size=0): polygon([(x0, y0), (x1, y0), (x1, y1), (x0, y1)], size) def label(x, y, material, size, group=0, circuit='', turns=0, angle=0): mi_addblocklabel(x, y) mi_selectlabel(x, y) mi_setblockprop(material, 0, size, circuit, angle, group, turns) mi_clearselected() def begin(axis=True, frequency=0, depth=20): global NODES, EDGES NODES, EDGES = {}, {} newdocument(0) mi_probdef(frequency, 'millimeters', 'axi' if axis else 'planar', 1e-10, depth, 25) mi_addmaterial('Воздух', 1, 1, 0, 0, 0) mi_addmaterial('Обмотка без потерь', 1, 1, 0, 0, 0) mi_addmaterial('Сталь линейная учебная', 1000, 1000, 0, 0, 0) mi_addmaterial('Сталь с учебной кривой насыщения', 1, 1, 0, 0, 0) for b, h in BH: mi_addbhpoint('Сталь с учебной кривой насыщения', b, h) mi_addboundprop('Внешняя граница A=0', 0, 0, 0, 0, 0, 0, 0, 0, 0) def exterior(axis=True, scale=1, mesh=4, cuts=()): bound = 'Внешняя граница A=0' radius, bottom, top = 120 * scale, -120 * scale, 120 * scale if axis: levels = sorted(set([bottom, top] + list(cuts))) for a, b in zip(levels[:-1], levels[1:]): edge((0, a), (0, b), 0, bound) for a, b in [((0, bottom), (radius, bottom)), ((radius, bottom), (radius, top)), ((radius, top), (0, top))]: edge(a, b, 0, bound) label(radius * .8, 0, 'Воздух', mesh) else: for a, b in [((-radius, bottom), (radius, bottom)), ((radius, bottom), (radius, top)), ((radius, top), (-radius, top)), ((-radius, top), (-radius, bottom))]: edge(a, b, 0, bound) label(radius * .8, 0, 'Воздух', mesh) def coil(r0, r1, z0, z1, name, turns, current, mesh, group=3): mi_addcircprop(name, current, 1) rect(r0, z0, r1, z1, mesh) label((r0 + r1) / 2, (z0 + z1) / 2, 'Обмотка без потерь', mesh, group, name, turns) def solve(tag, circuits=(), force_group=0, probe=(4, 0), **parameters): path = OUT + '/' + tag + '.fem' mi_setcomment('Учебная адаптация магнитных полей. Группы: 0 воздух, 1 сердечник, ' '2 исследуемое тело, 3 и 4 обмотки. Размеры мм, результаты СИ. ' + tag) mi_saveas(path) mi_close() opendocument(path) mi_analyze(1) mi_loadsolution() mo_clearblock() mo_groupselectblock() row = {'tag': tag, 'model': tag + '.fem', 'solution': tag + '.ans', 'elements': mo_numelements(), 'energy_J': real(mo_blockintegral(2)), 'coenergy_J': real(mo_blockintegral(17)), 'loss_W': real(mo_blockintegral(4))} mo_clearblock() if force_group: mo_groupselectblock(force_group) row['force_N'] = real(mo_blockintegral(19)) row['lorentz_N'] = real(mo_blockintegral(12)) mo_clearblock() row['circuits'] = [] for name in circuits: current, voltage, psi = mo_getcircuitproperties(name) row['circuits'].append({'name': name, 'current_A': parts(current), 'voltage_V': parts(voltage), 'psi_Wb': parts(psi)}) if probe: values = mo_getpointvalues(*probe) row.update(Bx_T=real(values[1]), By_T=real(values[2]), Hx_Am=real(values[5]), Hy_Am=real(values[6])) row['Mx_Am'] = row['Bx_T'] / MU0 - row['Hx_Am'] row['My_Am'] = row['By_T'] / MU0 - row['Hy_Am'] row.update(parameters) mo_close() mi_close() ROWS.append(row) with open(OUT + '/progress.json', 'w') as handle: json.dump({'solved': len(ROWS), 'last': tag}, handle) print('SOLVED', len(ROWS), tag, row['elements']) return row def psi(row, index=0): return complex(*row['circuits'][index]['psi_Wb']) def energy_check(row): work = sum(real(complex(*c['current_A']).conjugate() * complex(*c['psi_Wb'])) for c in row['circuits']) check('Энергия и коэнергия ' + row['tag'], rel(row['energy_J'] + row['coenergy_J'], work, 1e-9), .01) def finish(identifier, derived, settings): result = {'id': identifier, 'passed': all(c['passed'] for c in CHECKS), 'settings': settings, 'records': ROWS, 'derived': derived, 'checks': CHECKS} with open(OUT + '/result.json', 'w') as handle: json.dump(result, handle) print('RESULT', identifier, result['passed'], len(ROWS), len(CHECKS)) assert result['passed'], 'Не пройдены проверки; изучите result.json.' # F04: отдельная копия геометрии plunger-coenergy-force; уточняется интеграл работы. def plunger(position, mesh=.35, domain=1, suffix=''): begin() coil(10, 16, -10, 10, 'Обмотка возбуждения', 200, 2, mesh) rect(0, position - 8, 8, position + 8, mesh) label(4, position, 'Сталь линейная учебная', mesh, 2) rect(0, -35, 45, 50) label(40, 0, 'Воздух', .6) exterior(scale=domain, mesh=5, cuts=(-35, position - 8, position + 8, 50)) row = solve('F04-x%g-h%g-d%g%s' % (position, mesh, domain, suffix), ['Обмотка возбуждения'], 2, (4, position), position_mm=position, mesh_mm=mesh, domain=domain) row['L_H'] = real(psi(row)) / 2 energy_check(row) return row series = [plunger(x) for x in range(4, 21)] def integral(rows): return sum((b['position_mm'] - a['position_mm']) * .001 * (a['force_N'] + b['force_N']) / 2 for a, b in zip(rows[:-1], rows[1:])) work_fine, work_coarse = integral(series), integral(series[::2]) delta_coenergy = series[-1]['coenergy_J'] - series[0]['coenergy_J'] check('Шаг интегрирования работы 2/1 мм', rel(work_fine, work_coarse), .01) check('Работа и разность коэнергий при постоянном токе', rel(work_fine, delta_coenergy), .03) fine = plunger(12, .22) large = plunger(12, .22, 1.5) check('Сетка силы', rel(series[8]['force_N'], fine['force_N']), .03) check('Граница силы', rel(fine['force_N'], large['force_N']), .03) derivatives = [] for step in (.5, .25): low, high = plunger(12 - step, .16), plunger(12 + step, .16) derivatives.append((high['coenergy_J'] - low['coenergy_J']) / (.002 * step)) check('Шаг производной коэнергии', rel(*derivatives), .01) check('Сила по полю и коэнергии', rel(series[8]['force_N'], derivatives[-1]), .05) finish('plunger-work', {'work_J': work_fine, 'work_coarse_J': work_coarse, 'delta_coenergy_J': delta_coenergy, 'force_coenergy_N': derivatives[-1], 'series': [r['tag'] for r in series]}, {'N': 200, 'current_A': 2, 'coil_r_mm': [10, 16], 'coil_z_mm': [-10, 10], 'plunger_r_mm': 8, 'plunger_length_mm': 16, 'mu_r': 1000, 'position_definition': 'Координата центра плунжера z, а не длина воздушного зазора.', 'positive_force': 'В сторону увеличения z; на выбранном ходе сила отрицательна.'})