Взаимный нагрев кабелей и расстояние между ними

"""Самостоятельная лабораторная работа курса FEMM 5.
Запуск: встроенный MicroPython FEMM 5, команды FEMM доступны непосредственно.
Геометрические координаты и глубина — мм; результаты интегралов — единицы СИ.
Данные, допущения, аналитические ориентиры и ход исследования заданы в README.md.
Каждая серия записывается в новую calculated-XX; прежние результаты сохраняются.
Сохранённая модель закрывается и открывается до расчёта, включая русские имена.
"""

import os
import sys
import json
import math

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)
records = []
checks = []


def real(value):
    if isinstance(value, (tuple, list)):
        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 polygon(prefix, points):
    for x, y in points:
        globals()[prefix + "_addnode"](x, y)
    for i, (x, y) in enumerate(points):
        x2, y2 = points[(i + 1) % len(points)]
        globals()[prefix + "_addsegment"](x, y, x2, y2)


def geometry(prefix, job, mesh):
    if job.get("polygon"):
        polygon(prefix, job["polygon"])
    for x1, y1, x2, y2 in job.get("interfaces", []) + job.get("segments", []):
        globals()[prefix + "_addnode"](x1, y1)
        globals()[prefix + "_addnode"](x2, y2)
        globals()[prefix + "_addsegment"](x1, y1, x2, y2)
    for x1, y1, x2, y2, angle, name in job.get("arcs", []):
        globals()[prefix + "_addnode"](x1, y1)
        globals()[prefix + "_addnode"](x2, y2)
        globals()[prefix + "_addarc"](x1, y1, x2, y2, angle, 1)
        cx, cy = job.get("arc_center", (0, 0))
        mx = (x1 + x2) / 2 - cx
        my = (y1 + y2) / 2 - cy
        radius = math.sqrt((x1 - cx) ** 2 + (y1 - cy) ** 2)
        factor = radius / math.sqrt(mx * mx + my * my)
        globals()[prefix + "_selectarcsegment"](cx + mx * factor, cy + my * factor)
        if prefix == "ci":
            ci_setarcsegmentprop(1, "<None>", 0, 0, name or "<None>")
        else:
            hi_setarcsegmentprop(1, name or "<None>", 0, 0, "<None>")
        globals()[prefix + "_clearselected"]()
    for x, y, radius, name in job.get("circles", []):
        globals()[prefix + "_addnode"](x - radius, y)
        globals()[prefix + "_addnode"](x + radius, y)
        globals()[prefix + "_addarc"](x - radius, y, x + radius, y, 180, 1)
        globals()[prefix + "_addarc"](x + radius, y, x - radius, y, 180, 1)
        globals()[prefix + "_selectarcsegment"](x, y + radius)
        globals()[prefix + "_selectarcsegment"](x, y - radius)
        if prefix == "ci":
            ci_setarcsegmentprop(1, "<None>", 0, 0, name or "<None>")
        else:
            hi_setarcsegmentprop(1, name or "<None>", 0, 0, "<None>")
        globals()[prefix + "_clearselected"]()
    for x, y, name, group in job["blocks"]:
        globals()[prefix + "_addblocklabel"](x, y)
        globals()[prefix + "_selectlabel"](x, y)
        globals()[prefix + "_setblockprop"](
            name,
            0,
            mesh
            * job.get("block_mesh_multipliers", {}).get(
                str(group), job.get("mesh_multipliers", {}).get(name, 1)
            ),
            group,
        )
        globals()[prefix + "_clearselected"]()


def current(job, mesh, stem):
    newdocument(3)
    ci_probdef("millimeters", job.get("type", "planar"), 0, 1e-10, job["depth"], 30)
    for name, sx, sy in job["materials"]:
        ci_addmaterial(name, sx, sy, 1, 1, 0, 0)
    for name, voltage, amperes, kind in job["terminals"]:
        ci_addconductorprop(name, voltage, amperes, kind)
    geometry("ci", job, mesh)
    for x, y, name in job.get("electrode_segments", []):
        ci_selectsegment(x, y)
        ci_setsegmentprop("<None>", mesh, 0, 0, 0, name)
        ci_clearselected()
    path = directory + "/" + stem + ".fec"
    ci_saveas(path)
    ci_close()
    opendocument(path)
    ci_analyze(1)
    ci_loadsolution()
    co_smooth("off")
    terminals = {}
    terminal_power = 0
    balance = 0
    current_scale = 0
    for name, voltage, amperes, kind in job["terminals"]:
        u, i = co_getconductorproperties(name)
        terminals[name] = {"U_V": real(u), "I_A": real(i)}
        terminal_power += real(complex(u) * complex(i).conjugate())
        balance += complex(i)
        current_scale = max(current_scale, abs(complex(i)))
    powers = {}
    power = 0
    for x, y, name, group in job["blocks"]:
        if name == "<No Mesh>":
            continue
        co_groupselectblock(group)
        powers[str(group)] = real(co_blockintegral(0))
        power += powers[str(group)]
        co_clearblock()
    samples = []
    for x, y in job.get("samples", []):
        values = co_getpointvalues(x, y)
        samples.append(
            {
                "x_mm": x,
                "y_mm": y,
                "U_V": real(values[0]),
                "J_A_m2": math.sqrt(
                    abs(complex(values[1])) ** 2 + abs(complex(values[2])) ** 2
                ),
            }
        )
    row = {
        "model": stem + ".fec",
        "mesh_mm": mesh,
        "nodes": co_numnodes(),
        "elements": co_numelements(),
        "terminals": terminals,
        "P_W": power,
        "group_power_W": powers,
        "samples": samples,
        "current_balance_relative": abs(balance) / current_scale,
        "power_balance_relative": abs(power - terminal_power) / max(power, 1e-30),
    }
    check(stem + ": баланс токов", row["current_balance_relative"], 1e-6)
    check(stem + ": баланс мощности", row["power_balance_relative"], 1e-6)
    co_close()
    ci_close()
    return row


def heat(job, mesh, stem):
    newdocument(2)
    hi_probdef("millimeters", job.get("type", "planar"), 1e-10, job["depth"], 30)
    for name, kx, ky, q in job["materials"]:
        hi_addmaterial(name, kx, ky, q, 0)
    for boundary in job["boundaries"]:
        name, kind, temp, ambient, h = boundary[:5]
        beta = boundary[5] if len(boundary) > 5 else 0
        hi_addboundprop(name, kind, temp, 0, ambient, h, beta)
    geometry("hi", job, mesh)
    for x, y, name in job.get("boundary_segments", []):
        hi_selectsegment(x, y)
        hi_setsegmentprop(name, job.get("boundary_mesh", 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 x, y, name, group in job["blocks"]:
        ho_groupselectblock(group)
        averages[str(group)] = real(ho_blockintegral(0))
        volume = real(ho_blockintegral(2))
        q = [m[3] for m in job["materials"] if m[0] == name][0]
        generated += q * volume
        ho_clearblock()
    samples = []
    for x, y in job.get("samples", []):
        samples.append({"x_mm": x, "y_mm": y, "T_K": real(ho_getpointvalues(x, y)[0])})
    convection = 0
    for x1, y1, x2, y2, ambient, h in job.get("convection_edges", []):
        count = 200
        total = 0
        for n in range(count + 1):
            t = real(
                ho_getpointvalues(
                    x1 + (x2 - x1) * n / count, y1 + (y2 - y1) * n / count
                )[0]
            )
            total += (0.5 if n in (0, count) else 1) * (t - ambient)
        length = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) * 0.001
        convection += total / count * length * job["depth"] * 0.001 * h
    row = {
        "model": stem + ".feh",
        "mesh_mm": mesh,
        "nodes": ho_numnodes(),
        "elements": ho_numelements(),
        "samples": samples,
        "group_average_K": averages,
        "generated_W": generated,
        "convection_W": convection,
    }
    if job.get("convection_edges") and not job.get("skip_convection_balance"):
        check(
            stem + ": баланс конвективного теплоотвода",
            abs(convection - generated) / max(generated, 1e-30),
            0.005,
        )
    ho_close()
    hi_close()
    return row


JOB = {
    "polygon": [[-20000, -20000], [20000, -20000], [20000, 1000], [-20000, 1000]],
    "depth": 1000,
    "boundaries": [["Удалённая среда", 0, 293.15, 0, 0]],
    "boundary_segments": [
        [-20000, -1000, "Удалённая среда"],
        [20000, -1000, "Удалённая среда"],
        [0, -20000, "Удалённая среда"],
        [0, 1000, "Удалённая среда"],
    ],
    "boundary_mesh": 500,
    "mesh_multipliers": {"Грунт": 2000},
    "samples": [[0, 0]],
}
for level, mesh in enumerate((0.8, 0.4, 0.2)):
    readings = []
    for pitch in (150, 300, 600):
        job = json.loads(json.dumps(JOB))
        job["materials"] = [
            ["Медь", 380, 380, 20 / (math.pi * 0.005**2)],
            ["Изоляция", 0.2, 0.2, 0],
            ["Грунт", 1, 1, 0],
        ]
        job["circles"] = []
        job["interfaces"] = [
            [-1000, -1000, 1000, -1000],
            [-1000, -1000, -1000, 1000],
            [1000, -1000, 1000, 1000],
        ]
        job["boundary_segments"] += [
            [-10000, 1000, "Удалённая среда"],
            [10000, 1000, "Удалённая среда"],
        ]
        job["blocks"] = [[0, -500, "Грунт", 99], [5000, -5000, "Грунт", 98]]
        job["block_mesh_multipliers"] = {"99": 80, "98": 2000}
        for cable, x in enumerate((-pitch, 0, pitch)):
            job["circles"] += [
                [x, 0, 5, ""],
                [x, 0, 10, ""],
                [x, 0, 30, ""],
                [x, 0, 60, ""],
            ]
            job["blocks"] += [
                [x, 0, "Медь", 10 * cable + 1],
                [x + 7, 0, "Изоляция", 10 * cable + 2],
                [x + 20, 0, "Грунт", 10 * cable + 3],
                [x + 45, 0, "Грунт", 10 * cable + 4],
            ]
            job["block_mesh_multipliers"][str(10 * cable + 3)] = 5
            job["block_mesh_multipliers"][str(10 * cable + 4)] = 15
        result = heat(job, mesh, "mesh-%d-pitch-%d" % (level + 1, pitch))
        p = pitch * 0.001
        self_heat = (
            20 / (2 * math.pi) * math.log(2 / 0.01)
            + 20 / (2 * math.pi * 0.2) * math.log(10 / 5)
            + 20 / (4 * math.pi * 380)
        )
        neighbor = 2 * 20 / (2 * math.pi) * math.log(math.sqrt(p * p + 4) / p)
        result["pitch_mm"] = pitch
        result["T_center_K"] = result["samples"][0]["T_K"]
        result["analytic_halfspace_K"] = 293.15 + self_heat + neighbor
        check(
            "Изображения линейных источников, шаг %d мм" % pitch,
            abs(result["T_center_K"] - result["analytic_halfspace_K"]),
            1,
        )
        check("Полная мощность трёх кабелей", abs(result["generated_W"] - 60), 0.02)
        readings.append(result)
    records.append(
        {"mesh_mm": mesh, "readings": readings, "model": readings[-1]["model"]}
    )
check(
    "Сходимость температуры средней прокладки, К",
    abs(
        records[-1]["readings"][1]["T_center_K"]
        - records[-2]["readings"][1]["T_center_K"]
    ),
    0.15,
)
actual = (
    records[-1]["readings"][0]["T_center_K"] - records[-1]["readings"][-1]["T_center_K"]
)
reference = (
    records[-1]["readings"][0]["analytic_halfspace_K"]
    - records[-1]["readings"][-1]["analytic_halfspace_K"]
)
check(
    "Влияние расстояния против аналитической разности, К", abs(actual - reference), 0.15
)

result = {
    "runtime": "FEMM 5; embedded MicroPython",
    "records": records,
    "checks": checks,
    "passed": all(item["passed"] for item in checks),
}
with open(directory + "/result.json", "w") as output:
    json.dump(result, output)
print(json.dumps(result))
assert result["passed"]