Температура подземного кабеля

"""Температура трёх кабелей в грунте.
Среда: встроенный 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 = True
CENTERS = (0, 300, 600)
AMBIENT_K = 283.15
H = 10
SOURCES = (30.3 / 2596.72e-6, 0, 13.35 / 5925.04e-6, 0, 2.1 / 4146.9e-6, 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
    hi_addmaterial("Однородный грунт", 1, 1, 0, 0)
    # Расчётная область: нижний полукруг R=5.5 м с центром (0, 1 м).
    hi_addnode(-5500, 1000)
    hi_addnode(5500, 1000)
    hi_addarc(-5500, 1000, 5500, 1000, 180, 2)
    hi_addsegment(-5500, 1000, 5500, 1000)
    hi_selectsegment(0, 1000)
    hi_setsegmentprop("Конвекция в воздух", 100 / (2**level), 0, 0, 0, "<None>")
    hi_clearselected()
    hi_addblocklabel(0, -500)
    hi_selectlabel(0, -500)
    hi_setblockprop("Однородный грунт", 0, 200 / (2**level), 99)
    hi_clearselected()
    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,
    }
    # Численное интегрирование конвекции вдоль поверхности грунта.
    count = 880
    total = 0
    for n in range(count + 1):
        temperature = real(ho_getpointvalues(-5500 + 11000 * n / count, 1000)[0])
        total += (0.5 if n in (0, count) else 1) * (temperature - AMBIENT_K)
    convection = total * 11 / count * H
    row["convection_W"] = convection
    row["heat_balance_relative"] = abs(convection - generated) / generated
    row["T_ground_over_middle_K"] = real(ho_getpointvalues(300, 1000)[0])
    check("Баланс тепла, сетка %d" % (level + 1), row["heat_balance_relative"], 0.01)
    records.append(row)
    ho_close()
    hi_close()
check(
    "Сходимость температуры жилы, К",
    abs(records[-1]["T_center_K"][0] - records[-2]["T_center_K"][0]),
    0.25,
)
check(
    "Температура средней жилы против учебного ориентира 84 C, К",
    abs(records[-1]["T_center_K"][1] - 357.15),
    3,
)
# Общие потери трёх оболочек 2.1 Вт/м: по 0.7 Вт/м на кабель.
check(
    "Суммарная мощность заданных источников, Вт",
    abs(records[-1]["generated_W"] - (3 * (30.3 + 13.35) + 2.1)),
    0.1,
)

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"]