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

# -*- coding: utf-8 -*-
"""
MARS-12S10P SSSR 轴向磁拉力仿真 (Motor-CAD 前台, PyMotorCAD 驱动)
================================================================
参照 pss\\maxcalculator\\motorcad_export.py 的做法: 变量名候选列表逐个尝试,
探测结果如实记入 results JSON (失败也记, 不掩盖)。

流程: 载入原 .mot → 另存时间戳副本 → 开力计算开关 → 空载(I=0)求解 →
      负载(原电流)求解 → 探测轴向力输出(变量+波形) → 解析交叉校核 → JSON。

用法:
  python axial_force_run.py [--quit] [--skip-noload]
    --quit         完成后关闭 Motor-CAD (默认保持前台打开供人工检查)
    --skip-noload  只跑负载工况
"""
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

# 轴向力输出变量候选 (未实测, 探测式; Units_Force=kN 需留意单位)
FORCE_VARS = [
    "AxialForce", "Axial_Force", "ForceAxial", "Force_Axial",
    "NetAxialForce", "AFMAxialForce", "RotorAxialForce", "Rotor_Axial_Force",
    "StatorAxialForce", "AxialForceMean", "MeanAxialForce",
    "AxialForceAverage", "AxialForce_Load", "AxialForce_OC",
]
# 轴向力波形图候选 (get_magnetic_graph_point)
FORCE_GRAPHS = [
    "AxialForceVsAngle", "AxialForce", "ForceAxial", "Axial Force",
    "Fz", "ForceZ", "Force (Axial)",
]
# 气隙磁密波形候选 (解析校核用)
BG_GRAPHS = [
    "AirgapFluxDensity", "Airgap Flux Density", "AirgapFluxDensityOC",
    "BAirgap", "FluxDensityAirgap",
]


def probe_var(mc, names):
    """按候选名读变量, 返回 (名, 值); 全失败 (None, None)。"""
    for n in names:
        try:
            return n, mc.get_variable(n)
        except Exception:
            continue
    return None, None


def probe_graph(mc, names, npoints):
    """按候选名读波形 (逐点), 返回 {name, x, y} 或 None。"""
    for n in names:
        try:
            x0, y0 = mc.get_magnetic_graph_point(n, 0)
        except Exception:
            continue
        xs, ys = [x0], [y0]
        for i in range(1, npoints):
            try:
                x, y = mc.get_magnetic_graph_point(n, i)
            except Exception:
                break
            xs.append(x)
            ys.append(y)
        return {"name": n, "x": xs, "y": ys}
    return None


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 run_case(mc, tag, results, npoints):
    print("== 工况 [%s]: do_magnetic_calculation ..." % tag)
    t0 = time.time()
    mc.do_magnetic_calculation()
    dt = time.time() - t0
    print("   求解耗时 %.1f s" % dt)
    case = {"solve_seconds": dt}

    n, v = probe_var(mc, ["TorqueValueAveragePerCycle", "AverageTorque",
                          "MeanTorque", "ShaftTorque"])
    case["avg_torque"] = {"variable": n, "value": v}
    print("   平均转矩: %s = %s" % (n, v))

    n, v = probe_var(mc, FORCE_VARS)
    case["axial_force_var"] = {"variable": n, "value": v}
    print("   轴向力变量: %s = %s" % (n, v))

    g = probe_graph(mc, FORCE_GRAPHS, npoints)
    if g:
        case["axial_force_graph"] = {"name": g["name"], "stats": stats(g["y"]),
                                     "x": g["x"], "y": g["y"]}
        print("   轴向力波形 [%s]: %s" % (g["name"], stats(g["y"])))
    else:
        case["axial_force_graph"] = None
        print("   [警告] 轴向力波形候选全部失败")

    g = probe_graph(mc, BG_GRAPHS, npoints)
    if g:
        ys = g["y"]
        b2_mean = sum(b * b for b in ys) / len(ys)
        case["airgap_B_graph"] = {"name": g["name"], "stats": stats(ys),
                                  "B2_mean": b2_mean}
        print("   气隙磁密波形 [%s]: %s, mean(B^2)=%.4f"
              % (g["name"], stats(ys), b2_mean))
    else:
        case["airgap_B_graph"] = None
        print("   [提示] 气隙磁密波形候选失败, 解析校核转 GUI 人工读数")

    results["case_" + tag] = case
    return case


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

    from ansys.motorcad.core import MotorCAD
    print("启动 Motor-CAD (前台) ...")
    mc = MotorCAD()
    results = {"when": ts, "source_mot": os.path.basename(MOT_SRC)}
    try:
        mc.load_from_file(MOT_SRC)
        out_mot = os.path.join(OUT_DIR, "MARS_SSSR_axialforce_%s.mot" % ts)
        mc.save_to_file(out_mot)
        print("工作副本: %s" % out_mot)
        results["work_mot"] = out_mot

        # ---- 模型关键参数回读 (如实入档) ----
        params = {}
        for label, names in [
                ("Slot_Number", ["Slot_Number"]),
                ("Pole_Number", ["Pole_Number"]),
                ("Airgap_mm", ["Airgap"]),
                ("Magnet_Thickness_mm", ["Magnet_Thickness"]),
                ("ShaftSpeed_rpm", ["ShaftSpeed", "Shaft_Speed"]),
                ("PeakCurrent_A", ["PeakCurrent", "Peak_Current"]),
                ("PhaseAdvance_deg", ["PhaseAdvance", "Phase_Advance"]),
                ("Stator_Lam_Dia_mm", ["Stator_Lam_Dia"]),
                ("Stator_Bore_mm", ["Stator_Bore"]),
                ("AFM_D_Rotor_mm", ["AFM_D_Rotor"]),
                ("TorquePointsPerCycle", ["TorquePointsPerCycle"]),
                ("Units_Force", ["Units_Force"]),
        ]:
            n, v = probe_var(mc, names)
            params[label] = v
            print("  %s: %s = %s" % (label, n, v))
        results["params"] = params
        i_load = float(params["PeakCurrent_A"] or 0.0)
        npoints = int(params["TorquePointsPerCycle"] or 30) + 1

        # ---- 打开电磁力计算开关 ----
        for var in ["ElectromagneticForcesCalc_Load",
                    "ElectromagneticForcesCalc_OC"]:
            try:
                mc.set_variable(var, True)
                print("  [OK] %s = True" % var)
            except Exception as e:
                print("  [警告] %s 设置失败: %s" % (var, e))

        # ---- 工况 A: 空载 (I=0, 磁钢对定子铁芯的静态轴向吸力) ----
        if not skip_noload:
            name_i, _ = probe_var(mc, ["PeakCurrent", "Peak_Current"])
            mc.set_variable(name_i, 0.0)
            print("  %s -> 0 (空载)" % name_i)
            run_case(mc, "noload", results, npoints)
            mc.set_variable(name_i, i_load)
            print("  %s 恢复 %.3f A" % (name_i, i_load))

        # ---- 工况 B: 负载 (模型自带电流) ----
        run_case(mc, "load", results, npoints)

        # ---- 解析交叉校核: F ≈ A_gap/(2μ0) · mean(B²) ----
        try:
            d_out = float(params["AFM_D_Rotor_mm"]) * 1e-3
            d_in = float(params["Stator_Bore_mm"]) * 1e-3
            area = math.pi / 4.0 * (d_out ** 2 - d_in ** 2)
            results["analytic"] = {"area_m2": area,
                                   "note": "F=A/(2mu0)*mean(B^2), B 取仿真气隙磁密"}
            for tag in ("noload", "load"):
                case = results.get("case_" + tag) or {}
                bg = case.get("airgap_B_graph")
                if bg and bg.get("B2_mean"):
                    f_est = area / (2.0 * MU0) * bg["B2_mean"]
                    results["analytic"]["F_est_%s_N" % tag] = f_est
                    print("  解析估算 F_%s ≈ %.1f N (A=%.5f m², mean(B²)=%.4f)"
                          % (tag, f_est, area, bg["B2_mean"]))
        except Exception as e:
            results["analytic"] = {"failed": str(e)}

        mc.save_to_file(out_mot)
        res_path = os.path.join(OUT_DIR, "results_axialforce_%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)
        return 0
    finally:
        if quit_after:
            try:
                mc.quit()
            except Exception:
                pass
        else:
            print("[提示] Motor-CAD 保持前台打开供检查; 自动化场景加 --quit。")


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