"""FEMM 5 embedded MicroPython port of David Meeker's 10Aug2006 example. Run inside FEMM 5 with this folder as the working directory. FEMM commands are native globals. This script does not import an external femm/COM client. Currents/voltages are PEAK complex phasors, with exp(+j omega t) convention. The original approximate constant 2x2 impedance update is preserved. An additional FINAL solve ensures the solution and accepted currents agree. """ import math import json import sys import os import time if sys.implementation.name != "micropython": raise RuntimeError("Run this script in FEMM 5 embedded MicroPython, not CPython.") N1 = 260 N2 = 90 V1_RMS = 120.0 ZLOAD = complex(10.0, 0.0) IX_PEAK = 0.12 MAX_ITER = 100 # Original script used 0.1 V peak. Tighten the same residual norm for checking. KVL_TOL_PEAK_V = 0.001 SQRT2 = math.sqrt(2.0) V_SOURCE = complex(SQRT2 * V1_RMS, 0.0) CASES = globals().get("RUN_CASES", [ ("base", "mytransformer.fem"), ("refined", "mytransformer_refined.fem"), ]) def conj(z): return complex(z.real, -z.imag) def pair(z): return [float(z.real), float(z.imag)] def realvalue(value): # Real-valued block integrals can be floats rather than complex objects. return float(value.real) if isinstance(value, complex) else float(value) def phasor(z): return {"peak_re_im": pair(z), "rms_re_im": pair(z / SQRT2), "peak_magnitude": abs(z), "rms_magnitude": abs(z) / SQRT2, "phase_deg": math.atan2(z.imag, z.real) * 180.0 / math.pi} def solve2(a, b, c, d, u, w): determinant = a * d - b * c if abs(determinant) <= 1e-14 * max(abs(a * d), abs(b * c), 1.0): raise RuntimeError("Singular approximate transformer/load impedance.") return ((d * u - b * w) / determinant, (a * w - c * u) / determinant) def residual(v1, v2, i2): u1, u2 = V_SOURCE - v1, -v2 - ZLOAD * i2 return u1, u2, math.sqrt(abs(u1)**2 + abs(u2)**2) def dump_json(path, value): handle = open(path, "w") handle.write(json.dumps(value)) handle.write("\n") handle.close() def row(handle, stage, k, i1, i2, v1, v2, u1, u2, norm, di1, di2): values = [i1.real, i1.imag, i2.real, i2.imag, v1.real, v1.imag, v2.real, v2.imag, u1.real, u1.imag, u2.real, u2.imag, norm, norm / SQRT2, di1.real, di1.imag, di2.real, di2.imag] handle.write("%s,%d,%s\n" % (stage, k, ",".join(["%.17g" % x for x in values]))) handle.flush() def losses(group=None): mo_clearblock() if group is None: mo_groupselectblock() else: mo_groupselectblock(group) loss3 = realvalue(mo_blockintegral(3)) loss4 = realvalue(mo_blockintegral(4)) loss6 = realvalue(mo_blockintegral(6)) mo_clearblock() return {"lamination_or_hysteresis_integral_W": loss3, "resistive_integral_W": loss4, "total_W": loss6} def run_case(name, model): started = time.time() csv = open(name + "_iterations.csv", "w") csv.write("stage,iteration,i1_peak_re_A,i1_peak_im_A,i2_peak_re_A,i2_peak_im_A," "v1_peak_re_V,v1_peak_im_V,v2_peak_re_V,v2_peak_im_V," "kvl1_peak_re_V,kvl1_peak_im_V,kvl2_peak_re_V,kvl2_peak_im_V," "kvl_norm_peak_V,kvl_norm_rms_V,delta_i1_re_A,delta_i1_im_A," "delta_i2_re_A,delta_i2_im_A\n") # With one argument FEMM's open() opens a model; mode argument is file I/O. open(model) # Magnetic commands are registered only after opening a magnetic model. try: native_current_command = mi_modifycircprop except NameError: raise RuntimeError("The native magnetic FEMM global API is unavailable after opening the model.") mi_saveas(name + "_solved.fem") mi_modifycircprop("Первичная обмотка", 1, IX_PEAK) mi_modifycircprop("Вторичная обмотка", 1, 0) mi_analyze(1) mi_loadsolution() probe1 = mo_getcircuitproperties("Первичная обмотка") probe2 = mo_getcircuitproperties("Вторичная обмотка") z11 = probe1[1] / IX_PEAK z12 = probe2[1] / IX_PEAK z21 = z12 z22 = z11 * (float(N2) / N1)**2 a, b, c, d = z11, z12, z21, z22 + ZLOAD i1, i2 = solve2(a, b, c, d, V_SOURCE, 0j) stopped_on_residual = False for k in range(1, MAX_ITER + 1): mo_close() mi_modifycircprop("Первичная обмотка", 1, i1) mi_modifycircprop("Вторичная обмотка", 1, i2) mi_analyze(1) mi_loadsolution() r1 = mo_getcircuitproperties("Первичная обмотка") r2 = mo_getcircuitproperties("Вторичная обмотка") v1, v2 = r1[1], r2[1] u1, u2, norm = residual(v1, v2, i2) di1, di2 = solve2(a, b, c, d, u1, u2) row(csv, "iteration", k, i1, i2, v1, v2, u1, u2, norm, di1, di2) print("%s iteration %d: KVL peak norm %.9g V" % (name, k, norm)) # Preserve original update-before-stop order; final solve below fixes # its mismatch between the last currents and the last FEMM solution. i1 += di1 i2 += di2 if norm < KVL_TOL_PEAK_V: stopped_on_residual = True break mo_close() mi_modifycircprop("Первичная обмотка", 1, i1) mi_modifycircprop("Вторичная обмотка", 1, i2) mi_analyze(1) mi_loadsolution() r1 = mo_getcircuitproperties("Первичная обмотка") r2 = mo_getcircuitproperties("Вторичная обмотка") read_i1, v1, flux1 = r1 read_i2, v2, flux2 = r2 u1, u2, norm = residual(v1, v2, read_i2) row(csv, "final", k, read_i1, read_i2, v1, v2, u1, u2, norm, 0j, 0j) csv.close() info = mo_getprobleminfo() all_loss = losses() primary_loss = losses(1) secondary_loss = losses(2) core_loss = losses(3) sinput = 0.5 * V_SOURCE * conj(read_i1) # Both currents enter their winding. Load current leaves secondary. iload = -read_i2 vload = v2 sload_terminal = 0.5 * vload * conj(iload) sload_impedance = 0.5 * ZLOAD * abs(iload)**2 swindings = 0.5 * (v1 * conj(read_i1) + v2 * conj(read_i2)) power_residual = sinput.real - sload_impedance.real - all_loss["total_W"] circuit_loss_residual = swindings.real - all_loss["total_W"] published_i1 = complex(2.0446, -0.1424) published_i2 = complex(-5.8061, 0.0173) reference_errors = [abs(read_i1 - published_i1) / abs(published_i1), abs(read_i2 - published_i2) / abs(published_i2)] scale_p = max(abs(sinput.real), 1.0) checks = { "final_kvl_converged": stopped_on_residual and norm <= KVL_TOL_PEAK_V, "final_currents_match_prescription": max(abs(read_i1 - i1), abs(read_i2 - i2)) <= 1e-10 * max(abs(i1), abs(i2), 1.0), "load_terminal_matches_impedance": abs(sload_terminal - sload_impedance) <= 0.001 * max(abs(sload_impedance), 1.0), "positive_input_load_losses": sinput.real > 0 and sload_impedance.real > 0 and all_loss["total_W"] >= 0, "loss_integrals_sum": abs(all_loss["total_W"] - all_loss["lamination_or_hysteresis_integral_W"] - all_loss["resistive_integral_W"]) <= 1e-8 * scale_p, "input_load_volume_power_balance": abs(power_residual) <= max(0.001, 1e-4 * scale_p), "winding_volume_loss_balance": abs(circuit_loss_residual) <= max(0.001, 1e-4 * scale_p), "published_currents_within_3_percent": max(reference_errors) < 0.03, "frequency_depth_units_preserved": abs(info[1] - 60.0) < 1e-12 and abs(info[2] - 0.041275) < 1e-12 and abs(info[3] - 0.0254) < 1e-12, } result = { "case": name, "source_model": model, "solved_model": name + "_solved.fem", "interpreter": sys.implementation.name, "working_directory": os.getcwd(), "elapsed_seconds": time.time() - started, "iteration_count": k, "parameters": {"n1": N1, "n2": N2, "source_voltage_rms_V": V1_RMS, "source_voltage_peak_re_im_V": pair(V_SOURCE), "load_impedance_re_im_ohm": pair(ZLOAD), "probe_current_peak_A": IX_PEAK, "kvl_tolerance_peak_V": KVL_TOL_PEAK_V, "original_kvl_tolerance_peak_V": 0.1, "harmonic_convention": "peak phasors, exp(+j omega t), average S=0.5 V conj(I)"}, "problem_info_type_frequency_depth_m_lengthunit_m": list(info), "nodes": int(mo_numnodes()), "elements": int(mo_numelements()), "approximate_z_transformer_re_im_ohm": [[pair(z11), pair(z12)], [pair(z21), pair(z22)]], "probe_no_load_voltage_ratio": phasor(probe2[1] / probe1[1]), "primary_current_A": phasor(read_i1), "secondary_winding_current_A": phasor(read_i2), "load_current_A": phasor(iload), "primary_voltage_V": phasor(v1), "load_voltage_V": phasor(vload), "primary_flux_linkage_Wb_turn": phasor(flux1), "secondary_flux_linkage_Wb_turn": phasor(flux2), "final_kvl1_peak_re_im_V": pair(u1), "final_kvl2_peak_re_im_V": pair(u2), "final_kvl_norm_peak_V": norm, "final_kvl_norm_rms_V": norm / SQRT2, "input_complex_power_re_im_VA": pair(sinput), "load_terminal_complex_power_re_im_VA": pair(sload_terminal), "load_impedance_complex_power_re_im_VA": pair(sload_impedance), "windings_complex_power_re_im_VA": pair(swindings), "input_power_W": sinput.real, "load_power_W": sload_impedance.real, "efficiency": sload_impedance.real / sinput.real, "input_power_factor": sinput.real / abs(sinput), "volume_losses": {"all": all_loss, "primary_group_1": primary_loss, "secondary_group_2": secondary_loss, "core_group_3": core_loss}, "input_minus_load_minus_volume_losses_W": power_residual, "winding_real_power_minus_volume_losses_W": circuit_loss_residual, "ideal_turns_ratio_n2_over_n1": float(N2) / N1, "ideal_secondary_voltage_rms_V": V1_RMS * N2 / N1, "loaded_voltage_ratio_magnitude": abs(vload) / abs(v1), "voltage_droop_from_ideal_fraction": 1.0 - abs(vload) / (SQRT2 * V1_RMS * N2 / N1), "published_reference": {"source_url": "https://www.femm.info/doku/doku.php?id=MyTransformer", "currents_peak_re_im_A": [pair(published_i1), pair(published_i2)], "relative_complex_current_errors": reference_errors, "note": "Historical FEMM example, rounded published currents; not exact solver equivalence."}, "checks": checks, "all_checks_passed": all(checks.values()), "limitations": "Nonlinear fundamental time-harmonic 2D planar model. Phi_h=0: no nonzero hysteresis angle was specified. Integral 3 can contain lamination eddy-current loss; do not label it measured hysteresis loss. No thermal, end-winding, waveform/THD or hardware validation.", } dump_json(name + "_result.json", result) summary = open(name + "_result.txt", "w") summary.write("FEMM 5 native MicroPython transformer validation\n") summary.write("Case: %s; iterations: %d; nodes: %d; elements: %d\n" % (name, k, result["nodes"], result["elements"])) summary.write("Final KVL peak norm: %.12g V\n" % norm) summary.write("I1 peak: %.12g %+.12gj A; I2 winding peak: %.12g %+.12gj A\n" % (read_i1.real, read_i1.imag, read_i2.real, read_i2.imag)) summary.write("Uload RMS: %.12g V; Pin: %.12g W; Pload: %.12g W; FEMM losses: %.12g W\n" % (abs(vload) / SQRT2, sinput.real, sload_impedance.real, all_loss["total_W"])) summary.write("Pin-Pload-volume_losses: %.12g W\n" % power_residual) for key in sorted(checks): summary.write("%s: %s\n" % (key, checks[key])) summary.write(result["limitations"] + "\n") summary.close() mo_close() mi_close() return result results = [] try: for case_name, model_name in CASES: results.append(run_case(case_name, model_name)) combined = {"interpreter": sys.implementation.name, "cases": results} if len(results) == 2 and results[0]["case"] == "base" and results[1]["case"] == "refined": base, fine = results changes = {} for field in ["input_power_W", "load_power_W"]: changes[field] = abs(fine[field] - base[field]) / max(abs(fine[field]), 1e-12) for field in ["primary_current_A", "secondary_winding_current_A", "load_voltage_V"]: x, y = base[field]["peak_re_im"], fine[field]["peak_re_im"] zx, zy = complex(x[0], x[1]), complex(y[0], y[1]) changes[field] = abs(zy - zx) / max(abs(zy), 1e-12) lp, lf = base["volume_losses"]["all"]["total_W"], fine["volume_losses"]["all"]["total_W"] changes["volume_losses_W"] = abs(lf - lp) / max(abs(lf), 1e-12) comparison_checks = { "mesh_grew": fine["elements"] > base["elements"] and fine["nodes"] > base["nodes"], "currents_voltage_within_1_percent": max(changes["primary_current_A"], changes["secondary_winding_current_A"], changes["load_voltage_V"]) < 0.01, "input_load_power_within_2_percent": max(changes["input_power_W"], changes["load_power_W"]) < 0.02, "loss_power_within_5_percent": changes["volume_losses_W"] < 0.05, } combined["mesh_comparison"] = {"relative_changes_using_refined_denominator": changes, "element_count_ratio": float(fine["elements"]) / base["elements"], "checks": comparison_checks, "note": "One refinement repeat is a sensitivity check, not a formal asymptotic convergence study."} combined["all_checks_passed"] = all(x["all_checks_passed"] for x in results) and all(comparison_checks.values()) else: combined["all_checks_passed"] = all(x["all_checks_passed"] for x in results) dump_json("RESULT.json", combined) print("Transformer checks passed:", combined["all_checks_passed"]) if not combined["all_checks_passed"]: raise RuntimeError("Transformer validation checks failed; inspect result JSON and CSV.") except Exception as exc: dump_json("RUN_ERROR.json", {"error": str(exc), "completed_cases": [x["case"] for x in results], "interpreter": sys.implementation.name}) raise