pcb-afm-simulation-system/axial_mag_pull-master/axial_mag_pull/axial_force_final.py

# -*- coding: utf-8 -*-
"""
MARS SSSR 轴向磁拉力 — 正式计算
=================================
第5轮探测确认: AFM 力数据在 3D lumped 力图, 命名沿用径向机惯例, 其 "Fr"
(法向力) 在 AFM 2.5D 展开模型中即轴向力:
  Fr_Rotor_OL_Lumped / Fr_Rotor_OC_Lumped (转子, 负载/空载)
  Fr_Stator_OL_Lumped / Fr_Stator_OC_Lumped (定子, 反作用)
节点: 转子 10 (36°步, 首尾重复共11点), 定子 12 (30°步, 共13点); 单位 N。
两个径向切片 (sec1 r=28.25mm, sec2 r=34.75mm) 分别读, 节点求和+切片求和
得净轴向力; 按时间步扫描得波形。

校核: (a) 定子合力 ≈ -转子合力; (b) Σ(Ft×r) ≈ 电磁转矩图;
      (c) 解析 F ≈ A/(2μ0)·mean(B²)。

用法: python axial_force_final.py [--quit]
"""
import json
import math
import os
import sys
import time

BASE = os.path.dirname(os.path.abspath(__file__))
MOT_SRC = os.path.join(BASE, "MARS-12S10P_SSSR_D76-C150_V5.0-0819.mot")
OUT_DIR = os.path.join(BASE, "output_motorcad")

MU0 = 4e-7 * math.pi
SEC_RADII_MM = [28.25, 34.75]      # AFM_SectionCentreRadius_Array
MAX_TSTEPS = 64
MAX_NODES = 40


def stats(ys):
    if not ys:
        return None
    return {"mean": sum(ys) / len(ys), "min": min(ys), "max": max(ys),
            "pk2pk": max(ys) - min(ys), "n": len(ys)}


def read_nodes(mc, graph, sec, tstep):
    """读某时间步的全部节点 (x=角度, y=力N); 首尾重复点保留由调用方处理。"""
    xs, ys = [], []
    for i in range(MAX_NODES):
        try:
            x, y = mc.get_magnetic_3d_graph_point(graph, sec, i, tstep)
        except Exception:
            break
        xs.append(x)
        ys.append(y)
    return xs, ys


def net_force_series(mc, graph):
    """净力时间序列: 对两切片、去重节点求和; 返回 (series, meta)。"""
    series = []
    meta = {"sections": {}}
    for tstep in range(MAX_TSTEPS):
        total = 0.0
        got = False
        for sec in (1, 2):
            xs, ys = read_nodes(mc, graph, sec, tstep)
            if not ys:
                continue
            got = True
            # 首尾重复 (0°与360°同一节点) 则去掉末点
            n_unique = len(ys) - 1 if (len(xs) > 1 and
                                       abs(xs[-1] - xs[0] - 360.0) < 1e-6) \
                else len(ys)
            total += sum(ys[:n_unique])
            if tstep == 0:
                meta["sections"][sec] = {"n_points": len(ys),
                                         "n_unique": n_unique, "x": xs}
        if not got:
            break
        series.append(total)
    return series, meta


def torque_from_ft(mc, graph):
    """t=0 时刻 Σ(Ft×r) 粗校核 (Nm)。"""
    tq = 0.0
    for sec, r_mm in zip((1, 2), SEC_RADII_MM):
        xs, ys = read_nodes(mc, graph, sec, 0)
        if not ys:
            return None
        n_unique = len(ys) - 1 if (len(xs) > 1 and
                                   abs(xs[-1] - xs[0] - 360.0) < 1e-6) \
            else len(ys)
        tq += sum(ys[:n_unique]) * (r_mm * 1e-3)
    return tq


def read_2d(mc, graph, maxpts=64):
    xs, ys = [], []
    for i in range(maxpts):
        try:
            x, y = mc.get_magnetic_graph_point(graph, i)
        except Exception:
            break
        xs.append(x)
        ys.append(y)
    return xs, ys


def main(argv):
    quit_after = "--quit" in argv
    os.makedirs(OUT_DIR, exist_ok=True)
    ts = time.strftime("%m%d_%H%M%S")

    from ansys.motorcad.core import MotorCAD, set_motorcad_exe
    # Motor-CAD 2026R1 新装机器可能未注册 MOTORCAD_ACTIVEX (pymotorcad 0.8.x
    # 仍依赖它); 此时显式定位 exe, 可用环境变量 MOTORCAD_EXE 覆盖路径。
    if not os.environ.get("MOTORCAD_ACTIVEX"):
        exe = os.environ.get(
            "MOTORCAD_EXE",
            r"D:\Program Files\ANSYS Inc\v261\motorcad\MotorCAD.exe")
        if os.path.isfile(exe):
            set_motorcad_exe(exe)
    print("启动 Motor-CAD (前台) ...")
    mc = MotorCAD()
    # /SCRIPTING 模式在部分机器上窗口创建但不显示 (任务栏有图标点不开),
    # 强制可见; 已可见时无副作用 (2026-08-26 同事复现机实测该问题)
    try:
        mc.set_visible(True)
    except Exception as e:
        print("  [提示] set_visible 失败 (不影响计算): %s" % e)
    results = {"when": ts, "source_mot": os.path.basename(MOT_SRC),
               "convention_note": ("AFM 2.5D 展开模型中 Fr(法向)=轴向力; "
                                   "OL=负载(RMS 21A), OC=空载开路")}
    try:
        mc.load_from_file(MOT_SRC)
        out_mot = os.path.join(OUT_DIR, "MARS_SSSR_axialF_%s.mot" % ts)
        mc.save_to_file(out_mot)
        results["work_mot"] = out_mot

        for key, names in [("RMSCurrent_A", ["RMSCurrent"]),
                           ("ShaftSpeed_rpm", ["ShaftSpeed"]),
                           ("Airgap_mm", ["Airgap"]),
                           ("Stator_Lam_Dia_mm", ["Stator_Lam_Dia"]),
                           ("Stator_Bore_mm", ["Stator_Bore"])]:
            try:
                results[key] = mc.get_variable(names[0])
            except Exception:
                results[key] = None
        for var in ["ElectromagneticForcesCalc_Load",
                    "ElectromagneticForcesCalc_OC"]:
            mc.set_variable(var, True)

        print("求解 (负载点 RMS %sA, OC+OL 力同算) ..." % results["RMSCurrent_A"])
        t0 = time.time()
        mc.do_magnetic_calculation()
        results["solve_seconds"] = time.time() - t0
        print("  耗时 %.1f s" % results["solve_seconds"])

        # ---- 净轴向力: 转子/定子 x OL/OC ----
        forces = {}
        for graph in ["Fr_Rotor_OL_Lumped", "Fr_Stator_OL_Lumped",
                      "Fr_Rotor_OC_Lumped", "Fr_Stator_OC_Lumped"]:
            series, meta = net_force_series(mc, graph)
            if series:
                forces[graph] = {"series_N": series, "stats": stats(series),
                                 "meta": meta}
                print("  %s: %s" % (graph, stats(series)))
            else:
                forces[graph] = None
                print("  [如实] %s 无数据" % graph)
        results["axial_forces"] = forces

        # ---- 校核 a: 定转子合力反号 ----
        checks = {}
        for case in ("OL", "OC"):
            fr = forces.get("Fr_Rotor_%s_Lumped" % case)
            fs = forces.get("Fr_Stator_%s_Lumped" % case)
            if fr and fs:
                mr, ms = fr["stats"]["mean"], fs["stats"]["mean"]
                checks["action_reaction_%s" % case] = {
                    "rotor_mean_N": mr, "stator_mean_N": ms,
                    "imbalance_pct": abs(mr + ms) / max(abs(mr), 1e-9) * 100}
        # ---- 校核 b: Σ(Ft×r) vs 转矩 ----
        tq_ft = torque_from_ft(mc, "Ft_Rotor_OL_Lumped")
        _, tq_graph = None, None
        txs, tys = read_2d(mc, 17)          # id17 = 总转矩 (第4轮已辨认)
        tq_graph = stats(tys)["mean"] if tys else None
        checks["torque_crosscheck"] = {"sum_Ft_x_r_Nm_t0": tq_ft,
                                       "torque_graph_mean_Nm": tq_graph}
        # ---- 校核 c: 解析 F ≈ A/(2μ0)·mean(B²), B 取气隙磁密图 ----
        bxs, bys = read_2d(mc, "FluxDensityAirgap")
        if bys:
            b2 = sum(b * b for b in bys) / len(bys)
            d_out = float(results["Stator_Lam_Dia_mm"]) * 1e-3
            d_in = float(results["Stator_Bore_mm"]) * 1e-3
            area = math.pi / 4.0 * (d_out ** 2 - d_in ** 2)
            checks["analytic"] = {"mean_B2_T2": b2, "area_m2": area,
                                  "F_est_N": area / (2 * MU0) * b2}
        results["checks"] = checks
        print("校核: %s" % json.dumps(checks, ensure_ascii=False, indent=1))

        # ---- CSV 波形 ----
        csv_path = os.path.join(OUT_DIR, "axial_force_%s.csv" % ts)
        with open(csv_path, "w", encoding="utf-8") as f:
            f.write("tstep,Fr_Rotor_OL_N,Fr_Stator_OL_N,"
                    "Fr_Rotor_OC_N,Fr_Stator_OC_N\n")
            nmax = max(len(v["series_N"]) if v else 0
                       for v in forces.values())
            for i in range(nmax):
                row = [str(i)]
                for g in ["Fr_Rotor_OL_Lumped", "Fr_Stator_OL_Lumped",
                          "Fr_Rotor_OC_Lumped", "Fr_Stator_OC_Lumped"]:
                    v = forces.get(g)
                    row.append("%.4f" % v["series_N"][i]
                               if v and i < len(v["series_N"]) else "")
                f.write(",".join(row) + "\n")
        results["csv"] = csv_path

        res_path = os.path.join(OUT_DIR, "axialforce_final_%s.json" % ts)
        with open(res_path, "w", encoding="utf-8") as f:
            json.dump(results, f, ensure_ascii=False, indent=2)
        print("RESULTS: %s" % res_path)
        print("CSV: %s" % csv_path)

        # ---- 结论摘要 ----
        for case, label in (("OL", "负载(RMS 21A)"), ("OC", "空载")):
            v = forces.get("Fr_Rotor_%s_Lumped" % case)
            if v:
                s = v["stats"]
                print("[结论] %s 转子净轴向力: 均值 %.1f N, 纹波峰峰 %.2f N"
                      % (label, s["mean"], s["pk2pk"]))
        return 0
    finally:
        if quit_after:
            try:
                mc.quit()
            except Exception:
                pass
        else:
            print("[提示] Motor-CAD 保持前台打开供检查。")


if __name__ == "__main__":
    sys.exit(main(sys.argv[1:]))