pcb-afm-simulation-system/axial_mag_pull-master/axial_mag_pull/gui_app/solver.py

# -*- coding: utf-8 -*-
"""
轴向磁拉力求解核心 — 从 axial_force_final.py / axial_compare.py 提取,
封装为可被 GUI 调用的 Solver 类。支持单点仿真与气隙扫描两种模式。
所有日志通过 log_cb 回调输出到 GUI。
"""
import json
import math
import os
import time

MU0 = 4e-7 * math.pi
MAX_TSTEPS = 64
MAX_NODES = 40

# 全局 keepalive: keep_open=True 时保存 MotorCAD 对象引用, 防止 solver 被 GC
# 时连带销毁 MotorCAD → 关闭 Motor-CAD 进程 (GUI 模式下必须保持)
_KEEPALIVE_MC = 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 _read_nodes(mc, graph, sec, tstep):
    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 = {"sections": {}}
    for tstep in range(MAX_TSTEPS):
        total, got = 0.0, False
        for sec in (1, 2):
            xs, ys = _read_nodes(mc, graph, sec, tstep)
            if not ys:
                continue
            got = True
            nu = len(ys) - 1 if (len(xs) > 1 and
                                 abs(xs[-1] - xs[0] - 360.0) < 1e-6) else len(ys)
            total += sum(ys[:nu])
            if tstep == 0:
                meta["sections"][sec] = {"n_points": len(ys), "n_unique": nu}
        if not got:
            break
        series.append(total)
    return series, meta


def _torque_from_ft(mc, sec_radii_mm):
    tq = 0.0
    for sec, r_mm in zip((1, 2), sec_radii_mm):
        xs, ys = _read_nodes(mc, "Ft_Rotor_OL_Lumped", sec, 0)
        if not ys:
            return None
        nu = len(ys) - 1 if (len(xs) > 1 and
                             abs(xs[-1] - xs[0] - 360.0) < 1e-6) else len(ys)
        tq += sum(ys[:nu]) * (r_mm * 1e-3)
    return tq


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


class Solver:
    """轴向磁拉力求解器。参数通过构造函数传入, run_single() / run_sweep() 执行。"""

    def __init__(self, mot_path, out_dir, rms_current=21.0, speed_rpm=5000.0,
                 airgap_mm=1.0, magnet_temp_c=100.0, sec_radii_mm=None,
                 motorcad_exe=None, keep_open=True, log_cb=None):
        self.mot_path = mot_path
        self.out_dir = out_dir
        self.rms_current = rms_current
        self.speed_rpm = speed_rpm
        self.airgap_mm = airgap_mm
        self.magnet_temp_c = magnet_temp_c
        self.sec_radii_mm = sec_radii_mm or [28.25, 34.75]
        self.motorcad_exe = motorcad_exe
        self.keep_open = keep_open
        self.log = log_cb or (lambda msg: None)
        self.mc = None

    def _start_motorcad(self):
        from ansys.motorcad.core import MotorCAD, set_motorcad_exe
        if not os.environ.get("MOTORCAD_ACTIVEX") and self.motorcad_exe:
            if os.path.isfile(self.motorcad_exe):
                set_motorcad_exe(self.motorcad_exe)
                self.log("  未检测到 MOTORCAD_ACTIVEX, 已显式定位 exe")
        self.log("启动 Motor-CAD (前台) ...")
        self.mc = MotorCAD()
        try:
            self.mc.set_visible(True)
        except Exception as e:
            self.log("  [提示] set_visible 失败 (不影响计算): %s" % e)

    def _load_and_configure(self, tag):
        os.makedirs(self.out_dir, exist_ok=True)
        ts = time.strftime("%m%d_%H%M%S")
        self.mc.load_from_file(self.mot_path)
        out_mot = os.path.join(self.out_dir, "MARS_SSSR_%s_%s.mot" % (tag, ts))
        self.mc.save_to_file(out_mot)

        # 设置工况参数
        self.mc.set_variable("RMSCurrent", self.rms_current)
        self.mc.set_variable("ShaftSpeed", self.speed_rpm)
        self.mc.set_variable("Airgap", self.airgap_mm)
        self.mc.set_variable("Magnet_Temperature", self.magnet_temp_c)
        # 开力计算开关
        self.mc.set_variable("ElectromagneticForcesCalc_Load", True)
        self.mc.set_variable("ElectromagneticForcesCalc_OC", True)

        # 回读确认
        params = {}
        for var in ["RMSCurrent", "ShaftSpeed", "Airgap", "Magnet_Temperature",
                    "Stator_Lam_Dia", "Stator_Bore"]:
            try:
                params[var] = self.mc.get_variable(var)
            except Exception:
                params[var] = None
        self.log("  工况: RMS=%.1fA, %drpm, 气隙=%.2fmm, 磁钢%.0f°C"
                 % (params["RMSCurrent"], params["ShaftSpeed"],
                    params["Airgap"], params["Magnet_Temperature"]))
        return out_mot, params, ts

    def _solve_once(self, tag="axialF"):
        """执行一次求解, 返回 (results_dict, csv_path, json_path)"""
        out_mot, params, ts = self._load_and_configure(tag)
        self.log("求解中 (OC+OL 力同算) ...")
        t0 = time.time()
        self.mc.do_magnetic_calculation()
        dt = time.time() - t0
        self.log("  耗时 %.1f s" % dt)

        results = {"when": ts, "source_mot": os.path.basename(self.mot_path),
                   "work_mot": out_mot, "solve_seconds": dt,
                   "params": {"rms_current": params["RMSCurrent"],
                              "speed_rpm": params["ShaftSpeed"],
                              "airgap_mm": params["Airgap"],
                              "magnet_temp_c": params["Magnet_Temperature"]},
                   "convention_note": "AFM 2.5D 中 Fr(法向)=轴向力; 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(self.mc, graph)
            if series:
                forces[graph] = {"series_N": series, "stats": _stats(series),
                                 "meta": meta}
                s = _stats(series)
                self.log("  %s: 均值 %.1f N, 纹波 %.2f N"
                         % (graph, s["mean"], s["pk2pk"]))
            else:
                forces[graph] = None
                self.log("  [警告] %s 无数据" % graph)
        results["axial_forces"] = forces

        # 三判据
        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}
        tq_ft = _torque_from_ft(self.mc, self.sec_radii_mm)
        tys = _read_2d(self.mc, 17)
        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}
        bys = _read_2d(self.mc, "FluxDensityAirgap")
        if bys and params.get("Stator_Lam_Dia") and params.get("Stator_Bore"):
            b2 = sum(b * b for b in bys) / len(bys)
            d_out = float(params["Stator_Lam_Dia"]) * 1e-3
            d_in = float(params["Stator_Bore"]) * 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
        self.log("校核: %s" % json.dumps(checks, ensure_ascii=False))

        # CSV
        csv_path = os.path.join(self.out_dir, "axial_force_%s_%s.csv" % (tag, 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

        # JSON
        json_path = os.path.join(self.out_dir, "axialforce_%s_%s.json" % (tag, ts))
        with open(json_path, "w", encoding="utf-8") as f:
            json.dump(results, f, ensure_ascii=False, indent=2)

        # 结论摘要
        for case, label in (("OL", "负载(RMS %.0fA)" % self.rms_current),
                            ("OC", "空载")):
            v = forces.get("Fr_Rotor_%s_Lumped" % case)
            if v:
                s = v["stats"]
                self.log("[结论] %s 转子净轴向力: 均值 %.1f N, 纹波峰峰 %.2f N"
                         % (label, s["mean"], s["pk2pk"]))
        return results, csv_path, json_path

    def run_single(self):
        """单点仿真。返回 results dict。"""
        self._start_motorcad()
        try:
            results, csv_path, json_path = self._solve_once("single")
            self.log("结果文件: %s" % json_path)
            return results
        finally:
            if not self.keep_open:
                try:
                    self.mc.quit()
                except Exception:
                    pass
            else:
                global _KEEPALIVE_MC
                _KEEPALIVE_MC = self.mc  # 保持引用, 防止 GC 关闭 Motor-CAD
                self.log("[提示] Motor-CAD 保持前台打开供检查。")

    def run_sweep(self, gaps):
        """气隙扫描。gaps: 气隙列表(mm)。返回汇总 dict。"""
        self._start_motorcad()
        try:
            os.makedirs(self.out_dir, exist_ok=True)
            ts = time.strftime("%m%d_%H%M%S")
            self.mc.load_from_file(self.mot_path)
            out_mot = os.path.join(self.out_dir, "MARS_SSSR_sweep_%s.mot" % ts)
            self.mc.save_to_file(out_mot)

            self.mc.set_variable("RMSCurrent", self.rms_current)
            self.mc.set_variable("ShaftSpeed", self.speed_rpm)
            self.mc.set_variable("Magnet_Temperature", self.magnet_temp_c)
            self.mc.set_variable("ElectromagneticForcesCalc_Load", True)
            self.mc.set_variable("ElectromagneticForcesCalc_OC", True)

            sweep = {"when": ts, "magnet_temp_C": self.magnet_temp_c,
                     "rms_current_A": self.rms_current, "work_mot": out_mot,
                     "cases": []}
            for g in gaps:
                self.mc.set_variable("Airgap", g)
                back = self.mc.get_variable("Airgap")
                self.log("== 气隙 %.2f mm (回读 %s), 求解 ..." % (g, back))
                t0 = time.time()
                self.mc.do_magnetic_calculation()
                dt = time.time() - t0
                case = {"airgap_mm": back, "solve_seconds": dt}
                for graph, key in [("Fr_Rotor_OC_Lumped", "F_OC"),
                                   ("Fr_Rotor_OL_Lumped", "F_OL")]:
                    s, _ = _net_force_series(self.mc, graph)
                    case[key] = _stats(s)
                bys = _read_2d(self.mc, "FluxDensityAirgap")
                if bys:
                    case["B2_mean_T2"] = sum(b * b for b in bys) / len(bys)
                self.log("   F_OC=%.1f N, F_OL=%.1f N (耗时 %.0fs)"
                         % (case["F_OC"]["mean"], case["F_OL"]["mean"], dt))
                sweep["cases"].append(case)

            # 磁负刚度 (有限差分, OC 口径)
            cs = sweep["cases"]
            if len(cs) >= 2:
                f = [c["F_OC"]["mean"] for c in cs]
                g = [c["airgap_mm"] for c in cs]
                kneg = {}
                for i in range(len(cs) - 1):
                    kneg["seg_%.2f_%.2f" % (g[i], g[i + 1])] = \
                        -(f[i + 1] - f[i]) / (g[i + 1] - g[i])
                if len(cs) >= 3:
                    kneg["central"] = -(f[-1] - f[0]) / (g[-1] - g[0])
                sweep["kneg_N_per_mm"] = kneg
                self.log("磁负刚度: %s N/mm" % json.dumps(kneg, ensure_ascii=False))

            json_path = os.path.join(self.out_dir, "sweep_results_%s.json" % ts)
            with open(json_path, "w", encoding="utf-8") as f:
                json.dump(sweep, f, ensure_ascii=False, indent=2)
            self.log("扫描结果: %s" % json_path)
            return sweep
        finally:
            if not self.keep_open:
                try:
                    self.mc.quit()
                except Exception:
                    pass
            else:
                global _KEEPALIVE_MC
                _KEEPALIVE_MC = self.mc  # 保持引用, 防止 GC 关闭 Motor-CAD
                self.log("[提示] Motor-CAD 保持前台打开供检查。")