Тепловое сопротивление кабеля

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

import os
import sys
import math
import json

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, (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,
        }
    )


BURIED = False
CENTERS = (0,)
AMBIENT_K = 293.15
H = 5
SOURCES = (30 / (math.pi * 0.02875**2), 0, 0, 0, 0, 0)

RADII_MM = (28.75, 29.5, 52.5, 53, 57, 61)
NAMES = (
    "Медь жилы",
    "Экран жилы",
    "Маслобумажная изоляция",
    "Экран изоляции",
    "Свинцовая оболочка",
    "Защитный покров",
)
CONDUCTIVITIES = (380, 0.2, 0.2, 0.2, 35, 1 / 3.5)
# Теплопроводность оболочки принята равной 35 Вт/(м К).


def circle(x, radius, angle, boundary=""):
    hi_addnode(x - radius, 0)
    hi_addnode(x + radius, 0)
    hi_addarc(x - radius, 0, x + radius, 0, 180, angle)
    hi_addarc(x + radius, 0, x - radius, 0, 180, angle)
    hi_selectarcsegment(x, radius)
    hi_selectarcsegment(x, -radius)
    hi_setarcsegmentprop(angle, boundary or "<None>", 0, 0, "<None>")
    hi_clearselected()


for level, mesh in enumerate((3, 1.5, 0.75)):
    newdocument(2)
    hi_probdef("millimeters", "planar", 1e-10, 1000, 30)
    for number, name in enumerate(NAMES):
        k = CONDUCTIVITIES[number]
        hi_addmaterial(name, k, k, SOURCES[number], 0)
    hi_addboundprop("Конвекция в воздух", 2, 0, 0, AMBIENT_K, H, 0)
    angle = 4 / (2**level)
    for cable, x in enumerate(CENTERS):
        previous = 0
        for number, radius in enumerate(RADII_MM):
            boundary = "Конвекция в воздух" if not BURIED and number == 5 else ""
            circle(x, radius, angle, boundary)
            label_x = x + (previous + radius) / 2
            hi_addblocklabel(label_x, 0)
            hi_selectlabel(label_x, 0)
            hi_setblockprop(NAMES[number], 0, mesh, cable * 10 + number + 1)
            hi_clearselected()
            previous = radius

    model = directory + "/model-%02d.feh" % (level + 1)
    hi_saveas(model)
    hi_close()
    opendocument(model)
    hi_analyze(1)
    hi_loadsolution()
    ho_smooth("off")
    temperatures = [real(ho_getpointvalues(x, 0)[0]) for x in CENTERS]
    # Интеграл q_v dV использует фактическую площадь многоугольной сетки.
    generated = 0
    for cable, x in enumerate(CENTERS):
        for number, source in enumerate(SOURCES):
            if source:
                ho_groupselectblock(cable * 10 + number + 1)
                generated += source * real(ho_blockintegral(2))
                ho_clearblock()
    row = {
        "mesh_mm": mesh,
        "model": "model-%02d.feh" % (level + 1),
        "nodes": ho_numnodes(),
        "elements": ho_numelements(),
        "T_center_K": temperatures,
        "generated_W": generated,
    }
    interface = [real(ho_getpointvalues(r, 0)[0]) for r in RADII_MM]
    drops = [(temperatures[0] - interface[0]) / 30]
    drops += [(interface[i - 1] - interface[i]) / 30 for i in range(1, 6)]
    row["R_layers_K_m_W"] = drops
    row["T_interfaces_K"] = interface
    row["R_core_sheath_K_m_W"] = sum(drops[1:4])
    reference = []
    for i in range(1, 6):
        reference.append(
            math.log(RADII_MM[i] / RADII_MM[i - 1]) / (2 * math.pi * CONDUCTIVITIES[i])
        )
    row["R_layers_analytic_K_m_W"] = reference
    row["max_layer_relative_error"] = max(
        abs(drops[i + 1] / reference[i] - 1) for i in range(5)
    )
    check(
        "Радиальные сопротивления, сетка %d" % (level + 1),
        row["max_layer_relative_error"],
        0.03,
    )
    records.append(row)
    ho_close()
    hi_close()
check(
    "Сходимость температуры жилы, К",
    abs(records[-1]["T_center_K"][0] - records[-2]["T_center_K"][0]),
    0.25,
)
reference = 293.15 + 30 / (2 * math.pi * 0.061 * H)
reference += 30 / (4 * math.pi * 380)
reference += 30 * sum(records[-1]["R_layers_analytic_K_m_W"])
check(
    "Аналитическая температура центра, К",
    abs(records[-1]["T_center_K"][0] - reference),
    0.1,
)
check(
    "Сопротивления слоёв на окончательной сетке",
    records[-1]["max_layer_relative_error"],
    0.005,
)
check(
    "Мощность источника на окончательной сетке, Вт",
    abs(records[-1]["generated_W"] - 30),
    0.01,
)
check(
    "Сопротивление жила-оболочка против 0.488 К м/Вт",
    abs(records[-1]["R_core_sheath_K_m_W"] / 0.488 - 1),
    0.01,
)

result = {
    "runtime": "FEMM 5.0; 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"]