Индукционный нагрев: переходное тепловое поле и температурная обратная связь

"""Индукционный нагрев: нативная связь магнитного и теплового полей FEMM 5.

Запуск во встроенном MicroPython: команды mi/mo и hi/ho доступны глобально.
Геометрия — мм; температура — К; проводимость — См/м; потери — Вт;
объёмная теплоёмкость API hi — МДж/(м³ К). Ток катушки — амплитудный.
Мощность mo_blockintegral(4) уже усреднена по периоду; повторный множитель
1/2 при переносе тепла не применяется. Пространственная температура решается
нативно; в магнитные материалы возвращаются средние температуры ячеек.
Новый запуск создаёт calculated-NN рядом с программой. Прежние файлы
сохраняются. Модели закрываются и открываются повторно перед каждым расчётом.
"""

import os
import sys
import math
import json

CONFIG = {
    "frequency_Hz": 500,
    "coil_current_amplitude_A": 300,
    "coil_turns": 20,
    "sigma_0_S_m": 6e6,
    "alpha_sigma_K_inverse": 0.004,
    "mu_r_0": 20,
    "beta_mu_K_inverse": 0.005,
    "T_0_K": 293.15,
    "thermal_conductivity_W_mK": 25,
    "heat_capacity_MJ_m3K": 3.8,
    "convection_W_m2K": 1000,
    "end_time_s": 0.5,
    "iteration_temperature_tolerance_K": 0.0001,
    "iteration_limit": 12,
}
script = globals().get("__file__", "") or (sys.argv[0] if sys.argv else "")
base = os.path.dirname(os.path.abspath(script)) if script else os.getcwd()
index = 1
while "calculated-%02d" % index in os.listdir(base):
    index += 1
directory = base + "/calculated-%02d" % index
os.mkdir(directory)
os.chdir(directory)
records = []
checks = []


def real(value):
    if isinstance(value, (list, tuple)):
        return real(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 rectangle(prefix, x1, y1, x2, y2):
    points = [(x1, y1), (x2, y1), (x2, y2), (x1, y2)]
    for x, y in points:
        globals()[prefix + "_addnode"](x, y)
    for n, (x, y) in enumerate(points):
        xx, yy = points[(n + 1) % 4]
        globals()[prefix + "_addsegment"](x, y, xx, yy)


def cells(refined=False):
    radii = [0, 4, 7, 10]
    heights = [-2, 0, 2]
    if refined:
        radii = [0, 2, 4, 5.5, 7, 8.5, 10]
        heights = [-2, -1, 0, 1, 2]
    result = []
    for j in range(len(heights) - 1):
        for i in range(len(radii) - 1):
            a, b, c, d = radii[i], radii[i + 1], heights[j], heights[j + 1]
            result.append(
                {
                    "r1": a,
                    "r2": b,
                    "z1": c,
                    "z2": d,
                    "r": (a + b) / 2,
                    "z": (c + d) / 2,
                    "volume_m3": math.pi * (b * b - a * a) * (d - c) * 1e-9,
                }
            )
    return result


def grid(prefix, regions, mesh):
    # Одинаковая геометрия и ограничения сетки на всех шагах ANH-цепи.
    for cell in regions:
        rectangle(prefix, cell["r1"], cell["z1"], cell["r2"], cell["z2"])
    if prefix == "hi":
        for cell in regions:
            for x, y in (
                (cell["r"], cell["z1"]),
                (cell["r"], cell["z2"]),
                (cell["r1"], cell["z"]),
                (cell["r2"], cell["z"]),
            ):
                hi_selectsegment(x, y)
                hi_setsegmentprop("<None>", mesh, 0, 0, 0, "<None>")
                hi_clearselected()


def magnetic(regions, temperatures, mesh, stem):
    newdocument(0)
    smartmesh(0)
    mi_probdef(CONFIG["frequency_Hz"], "millimeters", "axi", 1e-9, 1, 30)
    mi_addmaterial("Воздух", 1, 1, 0, 0, 0)
    mi_addmaterial("Идеальная обмотка без собственных потерь", 1, 1, 0, 0, 0)
    mi_addboundprop("Удалённая граница A=0", 0, 0, 0, 0, 0, 0, 0, 0, 0)
    rectangle("mi", 0, -80, 80, 80)
    for x, y in ((0, 0), (80, 0), (40, -80), (40, 80)):
        mi_selectsegment(x, y)
        mi_setsegmentprop("Удалённая граница A=0", 0, 1, 0, 0)
        mi_clearselected()
    rectangle("mi", 0, -12, 24, 12)
    grid("mi", regions, mesh)
    conductivities = []
    permeabilities = []
    for n, cell in enumerate(regions):
        delta = temperatures[n] - CONFIG["T_0_K"]
        sigma = CONFIG["sigma_0_S_m"] / (1 + CONFIG["alpha_sigma_K_inverse"] * delta)
        mu = 1 + (CONFIG["mu_r_0"] - 1) / (1 + CONFIG["beta_mu_K_inverse"] * delta)
        assert sigma > 0 and mu > 1 and 250 < temperatures[n] < 600
        conductivities.append(sigma)
        permeabilities.append(mu)
        name = "Проводящая заготовка, ячейка %02d" % (n + 1)
        mi_addmaterial(name, mu, mu, 0, 0, sigma / 1e6)
        mi_addblocklabel(cell["r"], cell["z"])
        mi_selectlabel(cell["r"], cell["z"])
        mi_setblockprop(name, 0, mesh, "<None>", 0, 10 + n, 1)
        mi_clearselected()
    mi_addcircprop("Нагревательная катушка", CONFIG["coil_current_amplitude_A"], 1)
    rectangle("mi", 14, -4, 18, 4)
    for x, y, name, size, group, circuit, turns in (
        (
            16,
            0,
            "Идеальная обмотка без собственных потерь",
            2 * mesh,
            2,
            "Нагревательная катушка",
            CONFIG["coil_turns"],
        ),
        (22, 0, "Воздух", 2 * mesh, 0, "<None>", 1),
        (50, 0, "Воздух", 8 * mesh, 0, "<None>", 1),
    ):
        mi_addblocklabel(x, y)
        mi_selectlabel(x, y)
        mi_setblockprop(name, 0, size, circuit, 0, group, turns)
        mi_clearselected()
    path = directory + "/" + stem + ".fem"
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    mo_smooth("off")
    powers = []
    q = []
    for n, cell in enumerate(regions):
        mo_groupselectblock(10 + n)
        power = real(mo_blockintegral(4))
        volume = real(mo_blockintegral(10))
        mo_clearblock()
        check(
            stem + ": объём магнитной ячейки %d" % n,
            abs(volume / cell["volume_m3"] - 1),
            1e-7,
        )
        powers.append(power)
        q.append(power / volume)
    current, voltage, flux = mo_getcircuitproperties("Нагревательная катушка")
    terminal_power = 0.5 * (complex(voltage) * complex(current).conjugate()).real
    total = sum(powers)
    check(
        stem + ": амплитудная мощность выводов и потери",
        abs(total - terminal_power) / max(total, 1e-20),
        0.002,
    )
    independent = 0
    for i in range(1, mo_numelements() + 1):
        element = mo_getelement(i)
        group = int(element[6]) - 10
        if not 0 <= group < len(regions):
            continue
        x, y, area = element[3], element[4], element[5]
        value = mo_getpointvalues(x, y)
        density = (complex(value[7]) + complex(value[8])) * 1e6
        volume = area * 1e-6 * 2 * math.pi * x * 0.001
        independent += abs(density) ** 2 * volume / (2 * conductivities[group])
    check(
        stem + ": независимый интеграл J²/(2σ)",
        abs(independent - total) / max(total, 1e-20),
        0.03,
    )
    result = {
        "model": stem + ".fem",
        "elements": mo_numelements(),
        "power_W": total,
        "terminal_power_W": terminal_power,
        "independent_Joule_power_W": independent,
        "cell_q_W_m3": q,
        "cell_power_W": powers,
        "cell_material_T_K": temperatures[:],
        "cell_sigma_S_m": conductivities,
        "cell_mu_r": permeabilities,
    }
    mo_close()
    mi_close()
    return result


def thermal(regions, q, mesh, stem, previous="", dt=0):
    newdocument(2)
    smartmesh(0)
    hi_probdef("millimeters", "axi", 1e-10, 1, 30, previous, dt)
    t0 = CONFIG["T_0_K"]
    hi_addboundprop("Начальная температура", 0, t0, 0, 0, 0, 0)
    hi_addboundprop("Конвекция в среду", 2, 0, 0, t0, CONFIG["convection_W_m2K"], 0)
    grid("hi", regions, mesh)
    for n, cell in enumerate(regions):
        name = "Тепловая заготовка, ячейка %02d" % (n + 1)
        k = CONFIG["thermal_conductivity_W_mK"]
        hi_addmaterial(name, k, k, q[n], CONFIG["heat_capacity_MJ_m3K"])
        hi_addblocklabel(cell["r"], cell["z"])
        hi_selectlabel(cell["r"], cell["z"])
        hi_setblockprop(name, 0, mesh, 10 + n)
        hi_clearselected()
        for x, y in (
            (cell["r"], cell["z1"]),
            (cell["r"], cell["z2"]),
            (cell["r1"], cell["z"]),
            (cell["r2"], cell["z"]),
        ):
            if y not in (-2, 2) and x != 10:
                continue
            hi_selectsegment(x, y)
            hi_setsegmentprop(
                "Конвекция в среду" if previous else "Начальная температура",
                mesh,
                0,
                0,
                0,
                "<None>",
            )
            hi_clearselected()
    path = directory + "/" + stem + ".feh"
    hi_saveas(path)
    hi_close()
    opendocument(path)
    hi_analyze(1)
    hi_loadsolution()
    ho_smooth("off")
    averages = []
    generated = 0
    for n, cell in enumerate(regions):
        ho_groupselectblock(10 + n)
        averages.append(real(ho_blockintegral(0)))
        volume = real(ho_blockintegral(2))
        ho_clearblock()
        check(
            stem + ": объём тепловой ячейки %d" % n,
            abs(volume / cell["volume_m3"] - 1),
            1e-7,
        )
        generated += volume * q[n]
    heat_out = 0
    count = 160
    for a, b in (((0, -2), (10, -2)), ((0, 2), (10, 2)), ((10, -2), (10, 2))):
        total = 0
        for i in range(count + 1):
            fraction = i / count
            r = a[0] + (b[0] - a[0]) * fraction
            z = a[1] + (b[1] - a[1]) * fraction
            temperature = real(ho_getpointvalues(r, z)[0])
            total += (
                (0.5 if i in (0, count) else 1)
                * (temperature - t0)
                * 2
                * math.pi
                * r
                * 0.001
            )
        length = math.sqrt((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2) * 0.001
        heat_out += total / count * length * CONFIG["convection_W_m2K"]
    sampled = []
    for i in range(21):
        for j in range(9):
            r = (i + 0.001) / 20.002 * 10
            z = -2 + (j + 0.001) / 8.002 * 4
            sampled.append(real(ho_getpointvalues(r, z)[0]))
    average = sum(averages[n] * c["volume_m3"] for n, c in enumerate(regions)) / sum(
        c["volume_m3"] for c in regions
    )
    result = {
        "model": stem + ".feh",
        "previous_solution": previous,
        "elements": ho_numelements(),
        "cell_T_K": averages,
        "T_average_K": average,
        "T_sampled_max_K": max(sampled),
        "T_sampled_min_K": min(sampled),
        "heat_source_W": generated,
        "convection_W": heat_out,
    }
    ho_close()
    hi_close()
    return result


def run_series(mesh, dt, series, refined=False, frozen=False):
    regions = cells(refined)
    initial = thermal(regions, [0] * len(regions), mesh, "series-%d-initial" % series)
    previous = initial["model"].replace(".feh", ".anh")
    temperatures = [CONFIG["T_0_K"]] * len(regions)
    time = 0
    step = 0
    history = []
    max_balance = 0
    while time < CONFIG["end_time_s"] - 1e-10:
        step += 1
        actual_dt = min(dt, CONFIG["end_time_s"] - time)
        guess = temperatures[:]
        for iteration in range(1, CONFIG["iteration_limit"] + 1):
            stem = "series-%d-step-%02d-iteration-%02d" % (series, step, iteration)
            em = magnetic(
                regions,
                [CONFIG["T_0_K"]] * len(regions) if frozen else guess,
                mesh,
                stem + "-magnetic",
            )
            heat = thermal(
                regions, em["cell_q_W_m3"], mesh, stem + "-heat", previous, actual_dt
            )
            error = max(
                abs(heat["cell_T_K"][n] - guess[n]) for n in range(len(regions))
            )
            if frozen or error < CONFIG["iteration_temperature_tolerance_K"]:
                break
            guess = heat["cell_T_K"][:]
        else:
            raise AssertionError("Не сошлась связь магнитного и теплового полей")
        if not frozen:
            check(
                stem + ": сходимость температурной обратной связи, К",
                error,
                CONFIG["iteration_temperature_tolerance_K"],
            )
        check(
            stem + ": консервативный перенос потерь",
            abs(heat["heat_source_W"] / em["power_W"] - 1),
            1e-7,
        )
        storage = (
            CONFIG["heat_capacity_MJ_m3K"]
            * 1e6
            * sum(
                c["volume_m3"] * (heat["cell_T_K"][n] - temperatures[n]) / actual_dt
                for n, c in enumerate(regions)
            )
        )
        balance = abs(storage + heat["convection_W"] - heat["heat_source_W"]) / max(
            abs(storage), heat["heat_source_W"], 1e-20
        )
        check(stem + ": накопление энергии плюс конвекция", balance, 0.002)
        max_balance = max(max_balance, balance)
        time += actual_dt
        history.append(
            {
                "time_s": time,
                "T_average_K": heat["T_average_K"],
                "T_sampled_max_K": heat["T_sampled_max_K"],
                "T_sampled_min_K": heat["T_sampled_min_K"],
                "power_W": em["power_W"],
                "heat_out_W": heat["convection_W"],
                "storage_W": storage,
                "balance_relative": balance,
                "iterations": iteration,
                "feedback_residual_K": error,
                "magnetic_model": em["model"],
                "heat_model": heat["model"],
            }
        )
        temperatures = heat["cell_T_K"][:]
        penultimate = previous
        previous = heat["model"].replace(".feh", ".anh")
    check(
        "Серия %d: положительный пространственный нагрев" % series,
        max(0, CONFIG["T_0_K"] + 1 - heat["T_sampled_max_K"]),
        0,
    )
    return {
        "series": series,
        "mesh_mm": mesh,
        "dt_s": dt,
        "cells": len(regions),
        "frozen_materials": frozen,
        "time_s": time,
        "T_average_K": heat["T_average_K"],
        "T_sampled_max_K": heat["T_sampled_max_K"],
        "T_sampled_min_K": heat["T_sampled_min_K"],
        "final_power_W": em["power_W"],
        "magnetic_elements": em["elements"],
        "heat_elements": heat["elements"],
        "magnetic_model": em["model"],
        "heat_model": heat["model"],
        "previous_solution": penultimate,
        "maximum_balance_relative": max_balance,
        "cell_geometry": regions,
        "cell_T_K": temperatures,
        "cell_sigma_S_m": em["cell_sigma_S_m"],
        "cell_mu_r": em["cell_mu_r"],
        "cell_q_W_m3": em["cell_q_W_m3"],
        "history": history,
    }


for series, mesh in enumerate((0.8, 0.55, 0.35), 1):
    records.append(run_series(mesh, 0.125, series))
half = run_series(0.35, 0.0625, 4)
records.append(half)
refined = run_series(0.35, 0.0625, 5, refined=True)
records.append(refined)
frozen = run_series(0.35, 0.0625, 6, refined=True, frozen=True)
records.append(frozen)
fine = records[2]
scale = max(half["T_average_K"] - CONFIG["T_0_K"], 1)
check(
    "Чистое сгущение сетки: средняя температура",
    abs(fine["T_average_K"] - records[1]["T_average_K"]) / scale,
    0.015,
)
check(
    "Чистое сгущение сетки: мощность",
    abs(fine["final_power_W"] / records[1]["final_power_W"] - 1),
    0.025,
)
check(
    "Уточнение шага времени при неизменной сетке",
    abs(half["T_average_K"] - fine["T_average_K"]) / scale,
    0.025,
)
check(
    "Независимое уточнение переноса 6 → 24 ячейки",
    abs(refined["T_average_K"] - half["T_average_K"]) / scale,
    0.04,
)
feedback = abs(refined["final_power_W"] / frozen["final_power_W"] - 1)
check("Ненулевое влияние температурной обратной связи", max(0, 0.002 - feedback), 0)
check(
    "Пространственная неоднородность температуры разрешена",
    max(0, 0.5 - (refined["T_sampled_max_K"] - refined["T_sampled_min_K"])),
    0,
)
check(
    "Фактическое увеличение тепловой сетки",
    max(0, 1.5 - fine["heat_elements"] / records[0]["heat_elements"]),
    0,
)
check(
    "Фактическое увеличение магнитной сетки",
    max(0, 1.5 - fine["magnetic_elements"] / records[0]["magnetic_elements"]),
    0,
)
result = {
    "runtime": "FEMM 5; embedded MicroPython",
    "case_id": "native-induction-thermal-feedback",
    "records": records,
    "checks": checks,
    "passed": all(c["passed"] for c in checks),
    "plot": {
        "x": "series",
        "y": "T_average_K",
        "x_label": "Номер контрольной серии",
        "y_label": "Средняя температура, К",
    },
}
with open(directory + "/result.json", "w") as output:
    json.dump(result, output)
print(json.dumps(result))
assert result["passed"]