Источники тока, витки, знаки и связи между токами и напряжениями на выводах

"""Четыре витка: явная медь, MagnetWire и идеальный Litz при равной площади меди.
FEMM5/MicroPython; глубина 100 мм, I_peak=1 А в каждом витке.
Явные провода d=1 мм; MagnetWire: четыре таких витка в окне 2.5x2.5 мм;
Litz: четыре витка по 4 транспонированные жилы d=0.5 мм.
Разница Litz/явный одножильный — физическая замена провода, не ошибка FEMM.
"""

# Каталог сценария определяется явно: запуск не зависит от чужого CWD.
import os
import sys

os.chdir(os.path.dirname(os.path.abspath(sys.argv[0])))
exec(open("common.py").read())
folder = result_folder("results-winding-models")
rows = []
for frequency in [100, 10000]:
    for model in ["Explicit", "MagnetWire", "Litz"]:
        newdocument(0)
        mi_probdef(frequency, "millimeters", "planar", 1e-9, 100, 30)
        material("Air")
        mi_addcircprop("Winding", 1, 1)
        if model == "Explicit":
            material("Cu", sigma=58)
            for x in [-0.75, 0.75]:
                for y in [-0.75, 0.75]:
                    circle(x, y, 0.5)
                    label(x, y, "Cu", 0.04, "Winding", 1, 1)
        else:
            code = 3 if model == "MagnetWire" else 5
            strands = 1 if model == "MagnetWire" else 4
            diameter = 1 if model == "MagnetWire" else 0.5
            material("Cu", sigma=58, wire_type=code, strands=strands, diameter=diameter)
            mi_drawrectangle(-1.25, -1.25, 1.25, 1.25)
            label(0, 0, "Cu", 0.15, "Winding", 4, 1)
        label(3, 0, "Air", 1)
        mi_makeABC(7, 20, 0, 0, 0)
        solve(folder, "%s-%dHz" % (model, frequency))
        current, voltage, linkage = mo_getcircuitproperties("Winding")
        mo_groupselectblock(1)
        p4 = mo_blockintegral(4).real
        p6 = mo_blockintegral(6).real
        pin = (voltage * current.conjugate()).real / 2
        rows.append(
            {
                "model": model,
                "frequency_Hz": frequency,
                "P_ohmic_W": p4,
                "P_total_W": p6,
                "P_terminal_W": pin,
                "Z_ohm": pair(voltage / current),
                "Rdc_ohm": 4 * 0.1 / (58e6 * math.pi * 0.0005**2),
            }
        )
        close_magnetic()
        save_json(folder, rows)
assert all(row["P_total_W"] > 0 for row in rows)
print(rows)

"""Два официальных эталона, повторённых в FEMM5: точный момент и TEAM28.
Используются учебные torque-benchmark-course.fem и team28-course.fem с русскими именами свойств.
Диполь: m=1Ам²,B=1Тл; изменяется внутренний угол AGE без пересоздания сетки.
TEAM28: только установившаяся гармоническая сила при измеренной высоте11.3мм,
не полная переходная задача. Историческую силу автора не используем как ответ.
"""

# Каталог сценария определяется явно: запуск не зависит от чужого CWD.
import os
import sys

os.chdir(os.path.dirname(os.path.abspath(sys.argv[0])))
exec(open("common.py").read())
folder = result_folder("results-dipole-team28")
data = {"dipole": [], "team28": []}
for angle in [0, 30, 60, 90, 180, 270, 360]:
    opendocument(os.path.abspath("torque-benchmark-course.fem"))
    mi_modifyboundprop("Воздушный зазор AGE", 10, angle)
    solve(folder, "dipole-%d" % angle)
    torque = mo_gapintegral("Воздушный зазор AGE", 0)[0]
    exact = math.sin(angle * math.pi / 180)
    data["dipole"].append(
        {
            "angle_deg": angle,
            "torque_Nm": torque,
            "exact_Nm": exact,
            "absolute_error_Nm": abs(torque - exact),
            "elements": mo_numelements(),
        }
    )
    close_magnetic()
    save_json(folder, data)
for mesh_cm, boundary_name, radius_cm in [
    (0.125, "Asymptotic", 12),
    (0.0625, "Asymptotic", 12),
    (0.0625, "A0", 12),
    (0.0625, "A0", 24),
]:
    opendocument(os.path.abspath("team28-course.fem"))
    if boundary_name == "A0":
        # Меняется только физическое свойство внешней дуги.
        mi_modifyboundprop("Асимптотическая внешняя граница", 9, 0)
    if radius_cm != 12:
        # Увеличиваем только внешнюю дугу и два её осевых узла.
        # Катушки, диск, метки и заданная ближняя плотность сетки не меняются.
        mi_selectarcsegment(12, 0)
        mi_setgroup(99)
        mi_clearselected()
        mi_selectnode(0, -12)
        mi_selectnode(0, 12)
        mi_setgroup(99)
        mi_clearselected()
        mi_selectgroup(99)
        mi_scale(0, 0, radius_cm / 12, 4)
        mi_clearselected()
    # Размеры исходного файла в сантиметрах, ток20А — амплитуда.
    for r, z, name, turns in [
        (0.9, -0.1, "Воздух", 0),
        (4.2, -0.1, "Обмотка", 960),
        (8.8, 0.3, "Обмотка", -576),
        (3.89, 3.92, "Aluminum", 0),
    ]:
        # Имена исходного материала диска читаются непосредственно ниже.
        if name == "Aluminum":
            continue
        mi_selectlabel(r, z)
        mi_setblockprop(name, 0, mesh_cm, "Ток возбуждения" if turns else "<None>", 0, 0, turns)
        mi_clearselected()
    mi_selectlabel(3.89, 3.92)
    mi_setblockprop("Алюминиевый диск", 0, mesh_cm, "<None>", 0, 2, 0)
    mi_clearselected()
    solve(folder, "team28-%s-R%d-mesh-%g" % (boundary_name, radius_cm, mesh_cm))
    mo_selectblock(3.89, 3.92)
    data["team28"].append(
        {
            "mesh_cm": mesh_cm,
            "boundary": boundary_name,
            "outer_radius_cm": radius_cm,
            "Lorentz_Fz_N": mo_blockintegral(12).real,
            "WST_Fz_N": mo_blockintegral(19).real,
            "disc_power_W": mo_blockintegral(4).real,
            "weight_N": 0.107 * 9.80665,
            "elements": mo_numelements(),
        }
    )
    close_magnetic()
    save_json(folder, data)
assert max(x["absolute_error_Nm"] for x in data["dipole"]) < 0.002
print(data)

"""Проводящая труба: FE и независимое интегрирование точного радиального ОДУ.
FEMM5/MicroPython, без scipy. ri=12.7мм,ro=25.4мм,mu_r=1000,sigma=10МСм/м,
f=1Гц,I_peak=100А, глубина1м. Сначала труба I=0, затем заземлённая на бесконечности.
RK4 решает линейное комплексное ОДУ; 2000 и 4000 шагов проверяют интегрирование.
"""

# Каталог сценария определяется явно: запуск не зависит от чужого CWD.
import os
import sys

os.chdir(os.path.dirname(os.path.abspath(sys.argv[0])))
exec(open("common.py").read())
folder = result_folder("results-tube")
ri, ro = 0.0127, 0.0254
alpha2 = 1j * 2 * math.pi * 1000 * MU0 * 1e7


def integrate(h_initial, d_initial, count):
    step = (ro - ri) / count
    h, d = complex(h_initial), complex(d_initial)
    values = [[ri, h, d]]

    def slope(r, h, d):
        return d, (1 / r**2 + alpha2) * h - d / r

    for n in range(count):
        r = ri + n * step
        a, b = slope(r, h, d)
        c, e = slope(r + step / 2, h + a * step / 2, d + b * step / 2)
        f, g = slope(r + step / 2, h + c * step / 2, d + e * step / 2)
        k, l = slope(r + step, h + f * step, d + g * step)
        h += step * (a + 2 * c + 2 * f + k) / 6
        d += step * (b + 2 * e + 2 * g + l) / 6
        values.append([r + step, h, d])
    return values


def reference(count):
    a = integrate(100 / (2 * math.pi * ri), 0, count)
    b = integrate(0, 1, count)
    coefficient = (100 / (2 * math.pi * ro) - a[-1][1]) / b[-1][1]
    return [
        [v[0], v[1] + coefficient * w[1], v[2] + coefficient * w[2]]
        for v, w in zip(a, b)
    ]


ref = reference(4000)
ref_coarse = reference(2000)
ode_difference = max(abs(ref[2 * n][1] - ref_coarse[n][1]) for n in range(2001))
rows = []
for constrained, mesh in [(True, 0.8), (True, 0.4), (False, 0.4)]:
    newdocument(0)
    mi_probdef(1, "millimeters", "planar", 1e-10, 1000, 30)
    material("Air")
    material("Steel", 1000, 10)
    mi_addcircprop("Wire", 100, 1)
    mi_addcircprop("Tube", 0, 0)
    circle(0, 0, 2.54)
    circle(0, 0, 12.7)
    a0_circle(25.4)
    label(0, 0, "Air", mesh, "Wire", 1, 1)
    label(5, 0, "Air", mesh)
    label(20, 0, "Steel", mesh, "Tube" if constrained else "<None>", 1, 2)
    solve(folder, "Izero-%s-mesh-%g" % (constrained, mesh))
    samples = []
    numerator, denominator = 0, 0
    for index in range(100, 4000, 200):
        r, h, derivative = ref[index]
        exact_j = derivative + h / r
        p = mo_getpointvalues(r * 1000, 0)
        j = (p[7] + p[8]) * 1e6
        samples.append(
            {
                "r_m": r,
                "H_exact": pair(h),
                "H_FE": pair(p[6]),
                "J_exact": pair(exact_j),
                "J_FE": pair(j),
            }
        )
        numerator += abs(j - exact_j) ** 2
        denominator += abs(exact_j) ** 2
    mo_groupselectblock(2)
    rows.append(
        {
            "Izero": constrained,
            "mesh_mm": mesh,
            "tube_current_A": pair(mo_blockintegral(7)),
            "power_W": mo_blockintegral(4).real,
            "relative_RMS_J_error": math.sqrt(numerator / denominator),
            "ode_mesh_difference_A_per_m": ode_difference,
            "samples": samples,
        }
    )
    close_magnetic()
    save_json(folder, rows)
assert ode_difference < 1e-3
assert abs(complex(*rows[1]["tube_current_A"])) < 1e-5
assert rows[1]["relative_RMS_J_error"] < 0.08
print(rows)