# Машина с поверхностными постоянными магнитами # Шаг 11: Поле и численные результаты # Запускайте внутри FEMM 5 через File → Open Python Script. # Исходные учебные модели не изменяются; новый полевой расчёт не выполняется. import os try: _script = __file__.replace(chr(92), '/') _folder = _script.rsplit('/', 1)[0] if '/' in _script else os.getcwd() except NameError: _folder = os.getcwd() # Открываем готовый эталон и его унаследованное решение из одной папки. opendocument(_folder + '/' + "Поверхностные-магниты-09-Готовая-модель.fem") mi_loadsolution() mo_clearblock() mo_clearcontour() mo_zoomnatural() mo_showdensityplot(1, 0, 1.8, 0, 'bmag') mo_hidemesh() # Запрашиваем реальные результаты всех сохранённых цепей. print("Фаза А: I, U поля, потокосцепление", mo_getcircuitproperties("Фаза А")) print("Фаза В: I, U поля, потокосцепление", mo_getcircuitproperties("Фаза В")) print("Фаза С: I, U поля, потокосцепление", mo_getcircuitproperties("Фаза С")) # Энергия вычисляется по всем областям эталона. mo_groupselectblock() print('Вывод интеграла 2 (диагностика ошибки этой сборки), Дж:', mo_blockintegral(2)) mo_clearblock() # В этой радиальной модели FEMM 5.0 06Oct2026 игнорирует выражение направления # при вычислении H и энергии. Восстанавливаем направление для каждого элемента. # Используем B без сглаживания, площадь в единицах модели и глубину mo_getprobleminfo # в МЕТРАХ. После проверки сглаживание карты поля включается обратно. # Учебная редакция с русскими пояснениями; формулы и API исходника сохранены. """Энергия линейных изотропных постоянных магнитов в плоских моделях курса. Программу выполняют в Python FEMM либо передают словарь функций через api. Положения меток магнитов, Hc и возвратная проницаемость берутся из модели. Метод рассчитан на конкретные SPM/IPM-эталоны курса; нелинейная сталь, произвольное распределение намагниченности и 3D-эффекты сюда не входят.""" import math def integrate_linear_pm_energy(pm_labels, pm_mu, pm_hc, api=None): """Интегрировать несглаженные поля элементов и восстановить сглаживание on. pm_labels — список меток с x, y, direction и numeric_angle_deg. direction равно theta, theta+180 либо числовому направлению в градусах. Ближайшее угловое положение метки различает четыре разнесённых сектора. Материал магнитов в данных эталонах однозначно различается по pm_mu. pm_hc задаётся в А/м; pm_mu — относительная возвратная проницаемость. mo_getprobleminfo возвращает глубину в метрах, а mo_getelement — площадь в квадрате единиц модели: dV = area * scale**2 * depth_m, м³. У mi_getprobleminfo глубина, напротив, задаётся в единицах модели. Все возвращаемые интегральные энергии выражены в джоулях. При выходе сглаживание включается; прежнее состояние не сохраняется.""" if api is None: api = globals() # Проверяем плоскую магнитостатическую постановку до суммирования элементов. problem = api['mo_getprobleminfo']() if problem[0] != 0 or problem[1] != 0: raise ValueError('This helper requires a planar magnetostatic model') depth_m, scale = float(problem[2]), float(problem[3]) mu0 = 4 * math.pi * 1e-7 result = dict(raw_density_integral_J=0.0, physical_energy_J=0.0, coenergy_integral_J=0.0, correction_integral_J=0.0, pm_volume_m3=0.0, pm_element_count=0, Hparallel_min_A_m=1e300, Hparallel_max_A_m=-1e300) # Снимаем прежний выбор областей, чтобы интеграл относился ко всей модели. api['mo_clearblock']() api['mo_groupselectblock']() # Нативные интегралы служат независимой проверкой поэлементного суммирования. raw_energy = complex(api['mo_blockintegral'](2)).real raw_coenergy = complex(api['mo_blockintegral'](17)).real # Снимаем прежний выбор областей, чтобы интеграл относился ко всей модели. api['mo_clearblock']() # В центре элемента используем его собственное постоянное поле без сглаживания. api['mo_smooth']('off') try: for k in range(1, int(api['mo_numelements']()) + 1): element = api['mo_getelement'](k) x, y, area = map(float, element[3:6]) values = api['mo_getpointvalues'](x, y) bx, by = (complex(values[j]).real for j in (1, 2)) mux, muy = (complex(values[j]).real for j in (9, 10)) if mux <= 0 or abs(mux - muy) > 1e-8: raise ValueError('Only linear isotropic positive mu is supported') # Площадь переводится в м² отдельно от глубины, уже выраженной в метрах. volume = area * depth_m * scale ** 2 # Для линейной изотропной среды коэнергия имеет плотность B²/(2μ₀μr), Дж/м³. wc = (bx * bx + by * by) / (2 * mu0 * mux) w, correction = wc, 0.0 if abs(mux - pm_mu) < 1e-8: theta = math.atan2(y, x) * 180 / math.pi # Угловое расстояние приводим к интервалу −180…180° с учётом перехода через 0°. label = min(pm_labels, key=lambda lab: abs( (theta - math.atan2(lab['y'], lab['x']) * 180 / math.pi + 180) % 360 - 180)) direction = label['direction'] if direction == 'theta': alpha = theta elif direction == 'theta+180': alpha = theta + 180 else: alpha = float(direction) mx, my = math.cos(alpha * math.pi / 180), math.sin(alpha * math.pi / 180) # Числовой угол из метки нужен для сопоставления с нативным интегралом FEMM. raw_angle = label['numeric_angle_deg'] * math.pi / 180 rx, ry = math.cos(raw_angle), math.sin(raw_angle) # Остаточная индукция Br = μ₀μrHc; Hc — нормальная коэрцитивность. br = mu0 * mux * pm_hc # Энергия магнита вычисляется через B − Br·m; коэнергия и энергия здесь различаются. w = ((bx - br * mx) ** 2 + (by - br * my) ** 2) / (2 * mu0 * mux) # Поправка учитывает различие фактического и числового направлений намагничивания. correction = pm_hc * (bx * (rx - mx) + by * (ry - my)) # Проекция H вдоль намагничивания — локальная величина в А/м. hparallel = (bx * mx + by * my - br) / (mu0 * mux) result['Hparallel_min_A_m'] = min(result['Hparallel_min_A_m'], hparallel) result['Hparallel_max_A_m'] = max(result['Hparallel_max_A_m'], hparallel) result['pm_volume_m3'] += volume result['pm_element_count'] += 1 result['raw_density_integral_J'] += complex(values[4]).real * volume result['physical_energy_J'] += w * volume result['coenergy_integral_J'] += wc * volume result['correction_integral_J'] += correction * volume finally: api['mo_smooth']('on') result['raw_energy_J'] = raw_energy result['raw_coenergy_J'] = raw_coenergy result['corrected_from_raw_J'] = raw_energy + result['correction_integral_J'] agreements = dict( raw_density_vs_native_J=abs(result['raw_density_integral_J'] - raw_energy), coenergy_vs_native_J=abs(result['coenergy_integral_J'] - raw_coenergy), direct_vs_correction_J=abs(result['physical_energy_J'] - result['corrected_from_raw_J'])) # Допуск сравнения — абсолютный, в Дж, с масштабированием по энергии. tolerance_J = 1e-9 * max(1.0, abs(raw_energy), abs(raw_coenergy)) if max(agreements.values()) > tolerance_J: raise ValueError('Element integrals disagree: check units, materials and PM label identification') result['integral_agreement_abs_J'] = agreements result['integral_tolerance_J'] = tolerance_J result['integral_agreement_passed'] = True result['volume_formula'] = 'element_area * unit_scale**2 * postprocessor_depth_metres' result['method'] = 'unsmoothed planar element centroid sum; smoothing restored on' return result # Четыре метки именно приложенной модели; материал: mu_r=1.05, B_r=1.2 Тл. pm_labels = [{'x': 52.0, 'y': 0.0, 'direction': 'theta', 'numeric_angle_deg': 0.0}, {'x': 3.183976519938181e-15, 'y': 52.0, 'direction': 'theta+180', 'numeric_angle_deg': 0.0}, {'x': -52.0, 'y': 6.367953039876362e-15, 'direction': 'theta', 'numeric_angle_deg': 0.0}, {'x': -9.551929559814543e-15, 'y': -52.0, 'direction': 'theta+180', 'numeric_angle_deg': 0.0}] pm_energy = integrate_linear_pm_energy(pm_labels, 1.05, 1.2/(4*math.pi*1e-7*1.05)) print('Скорректированная энергия, Дж:', pm_energy['physical_energy_J']) print('Коэнергия, Дж:', pm_energy['coenergy_integral_J']) print('Проверка двух способов исправления, Дж:', pm_energy['physical_energy_J']-pm_energy['corrected_from_raw_J']) # Контроль физического поля на верхней оси: направление theta+180, m=(0,-1). # Контрольная точка использует сглаженное B, как исходная таблица поля. point = mo_getpointvalues(0, 52) hy_physical = (complex(point[2]).real + 1.2)/(4*math.pi*1e-7*1.05) print('Физическое Hy в точке (0;52) мм, А/м:', hy_physical) print('Проекция H вдоль магнитизации, А/м:', -hy_physical)