Индукционный нагрев трубы с зоной изменённой проводимости

"""FEMM 5: учебная модель, встроенный Python.
Команды mi_* и mo_* доступны непосредственно в FEMM; внешний пакет не нужен.
Повторный запуск создаёт отдельный calculated-NN и сохраняет предыдущие серии.
Геометрия в мм, токи и поля — комплексные амплитуды, результаты в СИ.
"""
import os
import sys
import math
import json

CONFIG = {'frequency_Hz': 5000, 'current_amplitude_A': 5,
          'turns': 50, 'conductivity_ratios': [1, 0.2],
          'outer_radius_mm': 600, 'check_outer_radius_mm': 800,
          'air_partition_radius_mm': 100}
OLD_FINAL = {
    '1': {'joint_loss_W': 2.4517742436608043, 'ohmic_loss_W': 6.875268455720551,
          'energy_J': 0.001215490293674114},
    '0.2': {'joint_loss_W': 1.4099030694851518, 'ohmic_loss_W': 8.103984233431396,
            'energy_J': 0.0012502526975481966}}
saved = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
base = os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd()
index = 1
while os.path.exists(base + '/calculated-%02d' % index):
    index += 1
directory = base + '/calculated-%02d' % index
os.mkdir(directory)
records = []
checks = []


def real(value):
    if isinstance(value, (list, tuple)):
        value = value[0]
    return complex(value).real


def check(name, value, limit):
    checks.append({'name': name, 'value': value, 'limit': limit,
                   'passed': math.isfinite(value) and value <= limit})


def material(name, conductivity=0):
    mi_addmaterial(name, 1, 1, 0, 0, conductivity / 1e6)


def label(r, z, name, mesh, group, circuit='<None>', turns=1):
    mi_addblocklabel(r, z)
    mi_selectlabel(r, z)
    mi_setblockprop(name, 0, mesh, circuit, 0, group, turns)
    mi_clearselected()


def rectangle(r1, z1, r2, z2):
    points = [(r1, z1), (r2, z1), (r2, z2), (r1, z2)]
    for r, z in points:
        mi_addnode(r, z)
    for i, point in enumerate(points):
        other = points[(i + 1) % 4]
        mi_addsegment(point[0], point[1], other[0], other[1])


def halfcircle(radius, angle, boundary=None):
    mi_addnode(0, -radius)
    mi_addnode(0, radius)
    mi_addarc(0, -radius, 0, radius, 180, angle)
    if boundary is not None:
        mi_selectarcsegment(radius, 0)
        mi_setarcsegmentprop(angle, boundary, 0, 0)
        mi_clearselected()


def solve(ratio, level, outer, boundary_check=False):
    scale = 2 ** level
    newdocument(0)
    mi_probdef(CONFIG['frequency_Hz'], 'millimeters', 'axi', 1e-10, 1000, 30, 0)
    mi_smartmesh(0)
    material('Воздух')
    material('Медь идеальной обмотки')
    material('Основной металл трубы', 6e6)
    material('Кольцевая зона трубы', 6e6 * ratio)
    mi_addboundprop('Удалённая внешняя граница', 0, 0, 0, 0, 0, 0, 0, 0, 0)
    halfcircle(outer, 4 / scale, 'Удалённая внешняя граница')
    mi_addsegment(0, -outer, 0, outer)
    # Только разделение сетки. На этой дуге нет граничного условия;
    # обе стороны имеют одинаковые mu_r=1 и sigma=0, поле непрерывно.
    halfcircle(CONFIG['air_partition_radius_mm'], 2 / scale)
    mi_addcircprop('Обмотка нагревателя', CONFIG['current_amplitude_A'], 1)
    rectangle(46, -10, 50, 10)
    label(48, 0, 'Медь идеальной обмотки', 1 / scale, 11,
          'Обмотка нагревателя', CONFIG['turns'])
    for lower, upper, name, group in [
            (-40, -4, 'Основной металл трубы', 5),
            (-4, 4, 'Кольцевая зона трубы', 6),
            (4, 40, 'Основной металл трубы', 7)]:
        rectangle(35, lower, 40, upper)
        label(37, (lower + upper) / 2, name, 0.5 / scale, group)
    label(0.5, 0, 'Воздух', 2 / scale, 3)
    label((outer + CONFIG['air_partition_radius_mm']) / 2, 0,
          'Воздух', 30 / scale, 1)
    suffix = '-R%d' % outer if boundary_check else ''
    name = 'model-%s-%02d%s.fem' % (str(ratio).replace('.', 'p'), level + 1, suffix)
    mi_saveas(directory + '/' + name)
    mi_close()
    opendocument(directory + '/' + name)
    mi_analyze(1)
    mi_loadsolution()
    mo_smooth('off')
    row = {'model': name, 'variant': ratio, 'level': level + 1,
           'outer_radius_mm': outer, 'boundary_check': boundary_check,
           'saved_reopened_solved': True,
           'nodes': mo_numnodes(), 'elements': mo_numelements()}
    mo_groupselectblock()
    row['ohmic_loss_W'] = real(mo_blockintegral(4))
    row['energy_J'] = real(mo_blockintegral(2))
    mo_clearblock()
    current, voltage, flux = mo_getcircuitproperties('Обмотка нагревателя')
    power = (complex(voltage) * complex(current).conjugate()).real / 2
    row['circuit_active_power_W'] = power
    row['power_balance_relative'] = abs(power / row['ohmic_loss_W'] - 1)
    check('Активная мощность цепи и потери: ' + name, row['power_balance_relative'], 0.002)

    expected_volume = {
        5: math.pi * (0.04 ** 2 - 0.035 ** 2) * 0.036,
        6: math.pi * (0.04 ** 2 - 0.035 ** 2) * 0.008,
        7: math.pi * (0.04 ** 2 - 0.035 ** 2) * 0.036,
        11: math.pi * (0.05 ** 2 - 0.046 ** 2) * 0.020}
    # Правая полуокружность в (r,z), вращаемая вокруг r=0, даёт целую сферу.
    expected_volume[3] = 4 * math.pi * 0.1 ** 3 / 3 - sum(expected_volume.values())
    expected_volume[1] = 4 * math.pi * ((outer * 0.001) ** 3 - 0.1 ** 3) / 3
    groups = [1, 3, 5, 6, 7, 11]
    counts = {str(group): 0 for group in groups}
    volumes = {str(group): 0 for group in groups}
    joule = 0
    # Независимое интегрирование по геометрии треугольников и J_phi.
    # mo_getpointvalues возвращает J в МА/м²; в СИ умножаем на 1e6.
    for element in range(1, int(mo_numelements()) + 1):
        values = mo_getelement(element)
        group = int(values[6])
        key = str(group)
        assert key in counts, 'Неожиданная группа элемента'
        counts[key] += 1
        r, z, area = values[3], values[4], values[5]
        volume = area * 1e-6 * 2 * math.pi * r * 0.001
        volumes[key] += volume
        if group in (5, 6, 7):
            point = mo_getpointvalues(r, z)
            density = (complex(point[7]) + complex(point[8])) * 1e6
            sigma = 6e6 * (ratio if group == 6 else 1)
            joule += abs(density) ** 2 * volume / (2 * sigma)
    row['element_counts_by_group'] = counts
    row['volumes_m3_by_group'] = volumes
    for group in groups:
        check('Использование каждой метки, группа%d: ' % group + name,
              0 if counts[str(group)] > 0 else 1, 0)
        error = abs(volumes[str(group)] / expected_volume[group] - 1)
        check('Объём области, группа%d: ' % group + name, error,
              0.003 / scale ** 2 if group in (1, 3) else 1e-8)
    row['element_Joule_integral_W'] = joule
    check('Независимый интеграл |J|²/(2σ): ' + name,
          abs(joule / row['ohmic_loss_W'] - 1), 0.02)
    losses = {}
    for group in (5, 6, 7):
        mo_groupselectblock(group)
        losses[str(group)] = real(mo_blockintegral(4))
        mo_clearblock()
    row['bulk_loss_W'] = losses['5'] + losses['7']
    row['joint_loss_W'] = losses['6']
    row['joint_loss_density_W_m3'] = losses['6'] / expected_volume[6]
    row['primary'] = losses['6']
    check('Сложение потерь трёх зон: ' + name,
          abs(sum(losses.values()) / row['ohmic_loss_W'] - 1), 1e-8)
    mo_close()
    mi_close()
    records.append(row)
    return row


for ratio in CONFIG['conductivity_ratios']:
    series = [solve(ratio, level, CONFIG['outer_radius_mm']) for level in range(3)]
    for key in ('joint_loss_W', 'ohmic_loss_W', 'energy_J'):
        check('Сходимость последних двух сеток, σз/σ=%s, %s' % (ratio, key),
              abs(series[2][key] / series[1][key] - 1), 0.012)
    control = solve(ratio, 2, CONFIG['check_outer_radius_mm'], True)
    for key in ('joint_loss_W', 'ohmic_loss_W', 'energy_J'):
        check('Удаление внешней границы600→800мм, σз/σ=%s, %s' % (ratio, key),
              abs(control[key] / series[2][key] - 1), 0.005)
    series[2]['correction_delta_relative_vs_v1'] = {
        key: series[2][key] / OLD_FINAL[str(ratio)] [key] - 1
        for key in ('joint_loss_W', 'ohmic_loss_W', 'energy_J')}

result = {'runtime': 'FEMM 5; embedded MicroPython',
          'case_id': 'pipe-joint-heating', 'model_version': 2,
          'records': records, 'checks': checks,
          'passed': all(item['passed'] for item in checks),
          'scope': 'Самостоятельная осесимметричная электромагнитная модель; вычисляется источник тепла, температура не рассчитывается.',
          'correction': 'R100мм разделяет только воздушную сетку. По одной метке в каждой области; физические материалы и граничное условие не изменены.'}
with open(directory + '/result.json', 'w') as stream:
    json.dump(result, stream)
print(json.dumps({'case_id': result['case_id'], 'model_version': 2,
                  'passed': result['passed'], 'models': len(records),
                  'checks': len(checks),
                  'failed': [item for item in checks if not item['passed']]}))
assert result['passed'], 'Численные проверки не пройдены; изучите result.json'