"""Учебное исследование магнитного поля во встроенном 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.' # F06: проводит сама сталь; отдельного кольца нет. N = 140 CURRENT = 2 * (20 - 10) * (35 - 10) / N COPPER_FILL = .5 COPPER_R = (2 * math.pi * .015 * N) / (58e6 * (250e-6 * COPPER_FILL / N)) def apparatus(frequency, gap, mesh=.35, domain=1, sigma=1, current=CURRENT): begin(frequency=frequency) mi_addmaterial('Массивная проводящая сталь', 100, 100, 0, 0, sigma) # Та же чашечная геометрия, что в bul-apparatus/pot-core. polygon([(0, 0), (30, 0), (30, 40), (22, 40), (22, 8), (8, 8), (8, 40), (0, 40)], mesh) rect(0, 40 + gap, 30, 48 + gap, mesh) label(4, 25, 'Массивная проводящая сталь', mesh, 1) label(15, 44 + gap, 'Массивная проводящая сталь', mesh, 2) coil(10, 20, 10, 35, 'Обмотка возбуждения', N, current, mesh) exterior(scale=domain, mesh=4, cuts=(0, 40, 40 + gap, 48 + gap)) row = solve('F06-f%g-g%g-h%g-d%g-s%g-i%g' % (frequency, gap, mesh, domain, sigma, current), ['Обмотка возбуждения'], 2, (4, 25), frequency_Hz=frequency, gap_mm=gap, mesh_mm=mesh, domain=domain, sigma_MSm=sigma, current_A=current) voltage = complex(*row['circuits'][0]['voltage_V']) impedance = voltage / current row['R_added_Ohm'] = impedance.real row['R_copper_Ohm'] = COPPER_R row['R_total_Ohm'] = impedance.real + COPPER_R row['L_equivalent_H'] = impedance.imag / (2 * math.pi * frequency) if frequency else real(psi(row)) / current row['loss_copper_W'] = COPPER_R * current ** 2 / (2 if frequency else 1) if frequency: check('Баланс мощности поля ' + row['tag'], rel(.5 * (voltage * current).real, row['loss_W'], 1e-8), .01) check('Пассивность ' + row['tag'], 0 if row['R_added_Ohm'] >= -1e-9 else 1, 0) return row series = [] for gap in (2, 6, 10): for frequency in (10, 50, 200, 500, 1000): row = apparatus(frequency, gap, .24) series.append(row['tag']) # Сетка и внешняя область в наиболее выраженном поверхностном эффекте. coarse = apparatus(1000, 2, .38) fine = apparatus(1000, 2, .17) large = apparatus(1000, 2, .17, 1.5) baseline = [r for r in ROWS if r['tag'] == series[4]][0] for key, limit in (('L_equivalent_H', .01), ('loss_W', .03), ('force_N', .03)): check('Сетка ' + key, rel(baseline[key], fine[key]), limit) check('Внешняя область ' + key, rel(fine[key], large[key]), limit) # Независимые контроли: отсутствие вихревых потерь и AC/DC нормировка силы. dc = apparatus(0, 6, .24, sigma=0) ac = apparatus(200, 6, .24, sigma=0) check('Средняя сила линейного AC равна половине DC при одинаковой амплитуде', rel(ac['force_N'], dc['force_N'] / 2), .01) check('Нулевая проводимость исключает вихревые потери', abs(ac['loss_W']), 1e-8) base = [r for r in ROWS if r['tag'] == series[7]][0] double = apparatus(200, 6, .24, current=2 * CURRENT) check('Квадратичность силы по току', rel(double['force_N'], 4 * base['force_N']), .01) check('Квадратичность потерь по току', rel(double['loss_W'], 4 * base['loss_W']), .01) finish('ac-solid-core', {'series': series}, {'N': N, 'current_peak_A': CURRENT, 'J_window_peak_A_mm2': 2, 'steel_mu_r': 100, 'steel_sigma_S_m': 1e6, 'copper_sigma_S_m': 58e6, 'copper_fill': COPPER_FILL, 'R_copper_Ohm': COPPER_R, 'minimum_skin_depth_mm': 1000 * math.sqrt(2 / (2 * math.pi * 1000 * MU0 * 100 * 1e6)), 'gap_definition': 'Два торцевых воздушных зазора между z=40 и нижней поверхностью дискового якоря z=40+gap.', 'force': 'Средняя осевая сила; притяжению соответствует отрицательный знак.', 'copper_model': 'Сосредоточенное DC-сопротивление отдельно от поля; торцевое удлинение и вытеснение в меди не учитываются.'})