From f42372ea68c0fe77c5d0e145031f34e28a736f9a Mon Sep 17 00:00:00 2001 From: smidl Date: Sun, 6 Sep 2026 11:07:47 +0100 Subject: [PATCH 1/2] Add the Python dq transform, its tests, and the Ansys/Python boundary rule Addresses the additive half of #25. The same Park/Clarke transform is written out by hand three times in `motors/synrm_5f_40s.py`'s output variables -- once for flux linkage, once for input current, once for terminal voltage -- off one pair of cos/sin tables. None of it needs the field solution; it is algebra on quantities the solve has already produced. This adds the Python side and the rule, and removes nothing: machine_design/transforms.py the transform once: to_dq / to_phases / electrical_angle, matching the Ansys convention exactly, including the sin(-h*(...)) sign that makes q the negative sine sum tests/test_transforms.py pins that convention offline (no Ansys) tests/test_ansys_transforms.py live cross-check on one solve, both subspaces excited CONTRIBUTING.md "What runs in Ansys, what runs in Python" Phase count is not hard-coded. `to_dq` reads it off the last axis of its input, so the same function serves the 3-phase anchor and the 5-phase one; `to_phases` takes it explicitly, having no input array to read it from. The offline tests are parametrised over (phase count, harmonic), and the m=3/h=3 triplen case -- zero-sequence, so not invertible -- is pinned as such so nobody "fixes" it. Verified on a live AEDT 2024.2 solve (bayes, 1 core): the Python transform reproduces Ansys's own dq output variables to 2.3e-15 (dq1) and 2.5e-16 (dq3) -- floating-point exact, not within a tolerance. Two things that verification turned up, both now recorded: * Ansys returns `theta_el` in DEGREES, not radians -- the `- pi` in its expression is evaluated as -180 in an angular context. The Ansys side is correct and internally consistent (theta_el is only ever used inside cos/sin as angle-minus-angle), but a Python consumer reading the raw number gets silent nonsense. `electrical_angle()` therefore takes an explicit `degrees=` flag and always returns radians. * `Flux_d1..q3` are defined as output variables but never extracted -- `solution_expressions` requests `Flux_e_*` (excitation flux) instead. So the flux transforms are computed by Ansys and discarded. The cross-check covers flux AND voltage deliberately. Flux stays in Ansys (it needs the solve), so that half is a permanent guard on the boundary. The voltage half is a one-shot equivalence proof: the Ansys-side voltage variables are slated to move to Python, which deletes both sides of the comparison, so it can only be run while they still exist. No behaviour change: three new files plus a CONTRIBUTING section. Co-Authored-By: Claude Opus 5 (1M context) --- CONTRIBUTING.md | 18 +++++ machine_design/transforms.py | 97 +++++++++++++++++++++++++ tests/test_ansys_transforms.py | 129 +++++++++++++++++++++++++++++++++ tests/test_transforms.py | 128 ++++++++++++++++++++++++++++++++ 4 files changed, 372 insertions(+) create mode 100644 machine_design/transforms.py create mode 100644 tests/test_ansys_transforms.py create mode 100644 tests/test_transforms.py diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index d809b9b..9d3a769 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -38,6 +38,24 @@ Each exports `Geometry` and `Computation`, so the module name carries the identi - **Derive, don't fork.** To study a variant (add magnets, change turns), subclass an existing anchor and override only what differs — do **not** edit an anchor's class in place, so old results stay reproducible. A genuinely new machine is a new anchor module. - **An anchor must declare all of the above**, rotor radii included. They are the radial window the barrier generators draw into; an anchor that does not set them cannot have a rotor built for it. +## What runs in Ansys, what runs in Python + +Ansys computes only what needs the field solution. Everything downstream of it is Python. + +| in Ansys (needs the solve) | in Python (algebra on solved quantities) | +|---|---| +| torque waveform, flux linkage, `InducedVoltage`, forces, core loss | Park/Clarke dq1/dq3 (`machine_design/transforms.py`), `Rs*I`, `Lew*dI/dt`, peak search over a waveform, terminal voltage at any speed | + +Why: + +- **One solve, many operating points.** Flux linkage at a given current is speed independent, so `V = Rs*I + j*w*psi` can be re-evaluated at any electrical speed in Python. An Ansys output variable is frozen at the solved frequency, so computing voltage there costs one solve *per* operating point -- three times the FEA for a three-point profile, N times for a drive cycle. +- **Post-processing scalars stay cheap.** `Lew`, winding turns, the bus voltage and `Rs` are post-processing quantities. Keeping them in Python means a stored pool of solves can be re-scored at a new value with no new FEA. +- **It can be tested.** Ansys expressions have no unit tests, parse differently across versions (#15), and mix units inside one argument (`2*pi*f*Time - 72deg`). The Python transforms are covered by `tests/test_transforms.py`. + +**One owner per quantity.** Where both sides *could* compute something, exactly one does, and the code says which. `Lew` in `motors/motor2.py` is the live example: it stays `0` in the model, and the end-winding voltage term belongs to the Python consumer. Setting it in both places double-counts it. + +`tests/test_ansys_transforms.py` (marked `ansys`) asserts that the Python transform and the Ansys output variables agree on a live solve, so this boundary is checked rather than merely described. + ## Pull requests - Shared code changes go through a branch and a PR, not a direct push to `main`. diff --git a/machine_design/transforms.py b/machine_design/transforms.py new file mode 100644 index 0000000..b59be76 --- /dev/null +++ b/machine_design/transforms.py @@ -0,0 +1,97 @@ +"""m-phase Park/Clarke transforms, in Python. + +This is the *single* Python implementation of the transform that Ansys applies in +`motors/synrm_5f_40s.py`'s output variables, where the same three lines are written out +once for flux linkage, once for input current and once for terminal voltage. The +convention here matches those expressions exactly, so a quantity transformed here +and the matching `*_d1`/`*_q1`/`*_d3`/`*_q3` output variable agree to solver +tolerance. `tests/test_ansys_transforms.py` asserts that on a live solve. + +Phase count and harmonic are both free. `to_dq` reads the phase count off the last +axis of its input, so the same function serves the 3-phase anchor and the 5-phase +one; `to_phases` takes it explicitly, since it has to synthesise that many values. +Caveat for the triplen case: in a 3-phase machine the 3rd harmonic is zero-sequence +(all phases in phase), so a dq3 read of it is not a constant vector -- that is +physics, not a defect of the transform. + +Convention (from `motors/synrm_5f_40s.py`, `set_output_vars`):: + + theta_el = RotSign * (Moving1.Position - InitPos) * PolePairs - pi + cos_ = cos( h * (theta_el - 2*pi*k/m)) + sin_ = sin(-h * (theta_el - 2*pi*k/m)) <-- note the minus + X_d = 2/m * sum_k X_k * cos_ + X_q = 2/m * sum_k X_k * sin_ + +(written out for m = 5 in that file; `m` here is the phase count.) + +**Units.** ``theta_el`` here is in RADIANS. Ansys's ``theta_el`` output variable comes +back in DEGREES -- the ``- pi`` in its expression is evaluated as -180 in an angular +context -- so convert with ``np.deg2rad`` before passing it in. Verified on a live +2024.2 solve: with that conversion the transform below reproduces Ansys's own +``Flux_d1/q1/d3/q3`` to 2e-15 (`tests/test_ansys_transforms.py`). + +The minus inside the sine means ``q`` is the *negative* of the usual sine sum. +That is a deliberate quirk of this model, not a typo; `tests/test_transforms.py` +pins it so it cannot drift. + +Why this lives in Python rather than as more Ansys output variables: see +"What runs in Ansys, what runs in Python" in `CONTRIBUTING.md`. In short, the +transform needs no field solution -- it is algebra on quantities the solve has +already produced -- and keeping it here makes it testable and lets one solve be +re-evaluated at any electrical speed. +""" + +import numpy as np + + +def electrical_angle(position, init_pos, pole_pairs, rot_sign=1.0, degrees=False): + """theta_el in RADIANS, matching Ansys's ``theta_el`` output variable. + + `position` and `init_pos` are mechanical angles -- radians by default, or degrees + with ``degrees=True``. Pass ``degrees=True`` for values read back from Ansys: + ``Moving1.Position`` is reported in degrees, so feeding it in as radians silently + produces a meaningless angle. The return value is always radians, which is what + `to_dq` and `to_phases` expect. + """ + position = np.asarray(position, dtype=float) + if degrees: + position, init_pos = np.deg2rad(position), np.deg2rad(init_pos) + return rot_sign * (position - init_pos) * pole_pairs - np.pi + + +def _phase_angles(theta_el, harmonic, n_phases): + """h * (theta_el - 2*pi*k/m) for k = 0..m-1, shaped (..., m).""" + theta_el = np.asarray(theta_el, dtype=float)[..., None] + k = np.arange(n_phases) + return harmonic * (theta_el - 2.0 * np.pi * k / n_phases) + + +def to_dq(x_phases, theta_el, harmonic): + """Per-phase quantity -> (d, q) for the given harmonic. + + `x_phases` carries the phases on the last axis, in winding order (A, B, C, ...); + the phase count is taken from that axis, so this works for any m. `theta_el` is + broadcast against the leading axes. Returns two arrays shaped like `theta_el`. + """ + x_phases = np.asarray(x_phases, dtype=float) + n_phases = x_phases.shape[-1] + if n_phases < 3: + raise ValueError(f"need at least 3 phases on the last axis, got {n_phases}") + phi = _phase_angles(theta_el, harmonic, n_phases) + scale = 2.0 / n_phases + d = scale * np.sum(x_phases * np.cos(phi), axis=-1) + q = scale * np.sum(x_phases * np.sin(-phi), axis=-1) + return d, q + + +def to_phases(d, q, theta_el, harmonic, n_phases): + """Inverse of `to_dq` for a single harmonic: (d, q) -> `n_phases` phase values. + + Exact inverse for a pure `harmonic` component; summing the results of several + harmonics reconstructs a multi-harmonic waveform, because the transforms for + distinct harmonics are orthogonal over m equally spaced phases (see + `tests/test_transforms.py`). `n_phases` is explicit here because, unlike `to_dq`, + there is no input array to read it from. + """ + phi = _phase_angles(theta_el, harmonic, n_phases) + return np.asarray(d, dtype=float)[..., None] * np.cos(phi) - np.asarray(q, dtype=float)[..., None] * np.sin(phi) diff --git a/tests/test_ansys_transforms.py b/tests/test_ansys_transforms.py new file mode 100644 index 0000000..9328355 --- /dev/null +++ b/tests/test_ansys_transforms.py @@ -0,0 +1,129 @@ +"""Live check that the Python dq transform agrees with the Ansys output variables. + +`machine_design/transforms.py` and `motors/synrm_5f_40s.py`'s `set_output_vars` implement the +same five-phase Park/Clarke transform twice -- once in Python, once as Ansys expressions. +This solves one rotor with both dq1 and dq3 excited, then asserts that transforming the +per-phase quantities in Python reproduces Ansys's own dq outputs. + +Two subjects, deliberately: + +* **flux linkage** -- `FluxLinkage(PhaseA..E)` vs `Flux_d1/q1/d3/q3`. Flux stays in Ansys + (it needs the field solution), so this comparison is permanent: it is the standing guard + on the boundary described in "What runs in Ansys, what runs in Python" (CONTRIBUTING.md). +* **terminal voltage** -- `V_A..V_E` vs `V_d1/q1/d3/q3`. This one is a *one-shot* + equivalence proof. The voltage transforms are slated to move to Python, which deletes + both sides of the comparison, so it can only be run while they still exist. Once the + Ansys-side voltage variables go, drop this half and keep the flux half. + +Marked `ansys`, so it is skipped by default. Run with: + + python -m pytest -m "" tests/test_ansys_transforms.py -v +""" + +import numpy as np +import pytest + +from machine_design.config import load_config +from machine_design.designs.design import Design +from machine_design.optimization.generators import HacklGenerator_OneLambda +from machine_design.transforms import to_dq +from motors.synrm_5f_40s import Computation, Geometry + +pytestmark = pytest.mark.ansys + +AEDT_VERSION = load_config()["aedt_version"] +R_STATOR_END = 0.7 +OFFSET = 0.35 +SEED = 42 +NUM_CORES = 1 + +# Both subspaces excited, so the h = 3 comparison is not a check against zero. +# Order is (Id1, Iq1, Id3, Iq3) -- compute() forwards these to set_variables(). +CURRENTS = (0.8, 0.9, 0.2, -0.15) + +PHASES = ["A", "B", "C", "D", "E"] +PHASE_FLUX = [f"FluxLinkage(Phase{p})" for p in PHASES] +PHASE_VOLTAGE = [f"V_{p}" for p in PHASES] + +# (label, per-phase expressions, ansys dq prefix) +SUBJECTS = [("flux", PHASE_FLUX, "Flux"), ("voltage", PHASE_VOLTAGE, "V")] + + +class ComputationWithPhaseQuantities(Computation): + """`Computation`, additionally extracting the per-phase quantities and theta_el. + + The shipped class extracts only the dq quantities, since that is all the optimizers + need. Subclassed rather than edited so the shared anchor is untouched. + """ + + def set_solution_expressions(self): + super().set_solution_expressions() + dq = [f"{prefix}_{axis}{h}" for prefix in ("Flux", "V") for h in (1, 3) for axis in "dq"] + extra = [*PHASE_FLUX, *PHASE_VOLTAGE, "theta_el", *dq] + self.solution_expressions = [*self.solution_expressions, *extra] + + +def _feasible_barriers(geometry): + dummy = Design(m2d=None, geometry=geometry, computation=Computation(geometry)) + generator = HacklGenerator_OneLambda(dummy, R_STATOR_END, offset=OFFSET) + + np.random.seed(SEED) + while True: + generator.set_parameters(generator.random_parameters()) + barriers = generator.split_barriers(generator.generate_barriers()) + if generator.feasible_barriers(barriers): + return barriers + + +@pytest.fixture(scope="module") +def solved(tmp_path_factory): + """One solve, shared by both subjects -- an FEA solve is far too costly to repeat.""" + geometry = Geometry() + computation = ComputationWithPhaseQuantities(geometry) + barriers = _feasible_barriers(geometry) + + design = None + try: + design = Design.create( + "TransformBoundary", + "Design01", + str(tmp_path_factory.mktemp("aedt") / "transform.aedt"), + geometry, + computation, + version=AEDT_VERSION, + non_graphical=True, + new_desktop=False, + close_on_exit=False, + ) + design.add_rotor() + for barrier in barriers: + design.add_rotor_barrier(barrier) + return design.compute(*CURRENTS, NUM_CORES=NUM_CORES) + finally: + if design is not None: + design.close_project() + + +@pytest.mark.parametrize(("label", "phase_exprs", "prefix"), SUBJECTS, ids=[s[0] for s in SUBJECTS]) +@pytest.mark.parametrize("harmonic", [1, 3]) +def test_python_transform_matches_ansys(solved, label, phase_exprs, prefix, harmonic): + assert solved is not None, "solve returned no results" + needed = [*phase_exprs, "theta_el", f"{prefix}_d{harmonic}", f"{prefix}_q{harmonic}"] + missing = [k for k in needed if k not in solved] + assert not missing, f"missing from the extracted solution: {missing}" + + # (n_timesteps, 5), phases in order A..E; to_dq reads the phase count off that axis. + per_phase = np.column_stack([np.asarray(solved[k], dtype=float) for k in phase_exprs]) + # Ansys returns theta_el in DEGREES (its "- pi" evaluates as -180 in an angular + # context); transforms.to_dq takes radians. + theta_el = np.deg2rad(np.asarray(solved["theta_el"], dtype=float)) + + d_py, q_py = to_dq(per_phase, theta_el, harmonic) + d_aedt = np.asarray(solved[f"{prefix}_d{harmonic}"], dtype=float) + q_aedt = np.asarray(solved[f"{prefix}_q{harmonic}"], dtype=float) + + # rtol handles solver/report round-off; atol keeps near-zero samples from failing on + # relative error alone, scaled to the signal rather than to absolute units. + atol = 1e-6 * max(np.max(np.abs(d_aedt)), np.max(np.abs(q_aedt)), 1.0) + np.testing.assert_allclose(d_py, d_aedt, rtol=1e-6, atol=atol, err_msg=f"{label} d{harmonic} mismatch") + np.testing.assert_allclose(q_py, q_aedt, rtol=1e-6, atol=atol, err_msg=f"{label} q{harmonic} mismatch") diff --git a/tests/test_transforms.py b/tests/test_transforms.py new file mode 100644 index 0000000..6fc8fd8 --- /dev/null +++ b/tests/test_transforms.py @@ -0,0 +1,128 @@ +"""Pins the m-phase dq convention used by `machine_design.transforms`. + +These need no Ansys. They exist because the same transform is written out by hand in +`motors/synrm_5f_40s.py`'s output variables, and a silent disagreement between the two +would corrupt every voltage and flux number downstream without failing loudly. The live +cross-check against Ansys is `tests/test_ansys_transforms.py`. + +Parametrised over phase count as well as harmonic, so the 3-phase anchor is covered by +the same transform as the 5-phase one. +""" + +import numpy as np +import pytest + +from machine_design.transforms import electrical_angle, to_dq, to_phases + +THETA = np.linspace(0.0, 2.0 * np.pi, 37) + +# (n_phases, harmonic). A harmonic that is a multiple of the phase count is +# zero-sequence and not invertible -- see test_triplen_harmonic_is_zero_sequence. +CASES = [(3, 1), (5, 1), (5, 3)] +CASE_IDS = [f"{m}ph-h{h}" for m, h in CASES] + + +@pytest.mark.parametrize(("n_phases", "harmonic"), CASES, ids=CASE_IDS) +@pytest.mark.parametrize(("d", "q"), [(1.0, 0.0), (0.0, 1.0), (0.7, -0.3), (-1.2, 2.5)]) +def test_round_trip_recovers_dq(n_phases, harmonic, d, q): + x = to_phases(d, q, THETA, harmonic, n_phases) + assert x.shape == (len(THETA), n_phases) + d_out, q_out = to_dq(x, THETA, harmonic) + np.testing.assert_allclose(d_out, d, atol=1e-12) + np.testing.assert_allclose(q_out, q, atol=1e-12) + + +@pytest.mark.parametrize(("h_signal", "h_read"), [(1, 3), (3, 1)]) +def test_harmonics_are_orthogonal(h_signal, h_read): + """A pure dq1 waveform must read as zero in dq3, and vice versa. + + This is what makes the two current subspaces independent design variables; if it + broke, dq3 injection would leak into the dq1 measurement. + """ + x = to_phases(1.3, -0.8, THETA, h_signal, 5) + d, q = to_dq(x, THETA, h_read) + np.testing.assert_allclose(d, 0.0, atol=1e-12) + np.testing.assert_allclose(q, 0.0, atol=1e-12) + + +def test_superposed_harmonics_separate(): + """dq1 and dq3 components ride on the same phase waveform without interfering.""" + x = to_phases(2.0, 0.5, THETA, 1, 5) + to_phases(-0.4, 0.9, THETA, 3, 5) + d1, q1 = to_dq(x, THETA, 1) + d3, q3 = to_dq(x, THETA, 3) + np.testing.assert_allclose(d1, 2.0, atol=1e-12) + np.testing.assert_allclose(q1, 0.5, atol=1e-12) + np.testing.assert_allclose(d3, -0.4, atol=1e-12) + np.testing.assert_allclose(q3, 0.9, atol=1e-12) + + +def test_triplen_harmonic_is_zero_sequence_in_three_phase(): + """m = 3, h = 3 is deliberately NOT invertible, and that is physics. + + 3*(2*pi*k/3) is a whole number of turns, so every phase sees the same angle: the + triplen harmonic is zero-sequence. Pinned so nobody "fixes" the transform to make + this round-trip. + """ + x = to_phases(1.0, 0.0, THETA, 3, 3) + np.testing.assert_allclose(x[:, 0], x[:, 1], atol=1e-12) + np.testing.assert_allclose(x[:, 0], x[:, 2], atol=1e-12) + + d, _ = to_dq(x, THETA, 3) + assert not np.allclose(d, 1.0, atol=1e-6), "expected the triplen round-trip to fail" + + +@pytest.mark.parametrize("n_phases", [3, 5]) +def test_amplitude_scaling_is_two_over_m(n_phases): + """The 2/m factor: a unit-amplitude d-axis waveform must give d = 1.""" + d, _ = to_dq(to_phases(1.0, 0.0, THETA, 1, n_phases), THETA, 1) + np.testing.assert_allclose(d, 1.0, atol=1e-12) + + +@pytest.mark.parametrize("n_phases", [3, 5]) +def test_q_axis_sign_convention(n_phases): + """`q` is the NEGATIVE sine sum, matching sin(-h*(theta - 2*pi*k/m)) in Ansys. + + Built from explicit cos/sin phase quantities rather than via `to_phases`, so this + fails if the sign in `to_dq` is ever "corrected" to the more usual convention. + """ + theta = 0.37 + phi = theta - 2.0 * np.pi * np.arange(n_phases) / n_phases + + # A pure +cos(phi) phase quantity is pure d-axis, unit amplitude. + d, q = to_dq(np.cos(phi), theta, 1) + assert d == pytest.approx(1.0, abs=1e-12) + assert q == pytest.approx(0.0, abs=1e-12) + + # A pure +sin(phi) phase quantity gives q = -1 under this convention, not +1. + d, q = to_dq(np.sin(phi), theta, 1) + assert d == pytest.approx(0.0, abs=1e-12) + assert q == pytest.approx(-1.0, abs=1e-12) + + +def test_electrical_angle_matches_ansys_expression(): + """theta_el = RotSign*(Position - InitPos)*PolePairs - pi, in radians.""" + position = np.array([0.0, 0.25, 0.5]) + init_pos = np.deg2rad(-45.0) + got = electrical_angle(position, init_pos, pole_pairs=2, rot_sign=1.0) + expected = 1.0 * (position - init_pos) * 2 - np.pi + np.testing.assert_allclose(got, expected, atol=1e-12) + + +def test_electrical_angle_accepts_degrees_like_ansys_reports_them(): + """Ansys reports `Moving1.Position` in degrees; `degrees=True` must match radians. + + Guards the trap this API invites: the natural input comes back from Ansys in + degrees, and passing it as radians fails silently rather than loudly. + """ + position_deg, init_deg = np.array([315.0, 316.8, 318.6]), -45.0 + from_deg = electrical_angle(position_deg, init_deg, pole_pairs=2, degrees=True) + from_rad = electrical_angle(np.deg2rad(position_deg), np.deg2rad(init_deg), pole_pairs=2) + np.testing.assert_allclose(from_deg, from_rad, atol=1e-12) + + # The observed Ansys value: (315 - -45) * 2 - 180 degrees = 540 deg = 3*pi rad. + assert from_deg[0] == pytest.approx(3.0 * np.pi, abs=1e-12) + + +def test_rejects_too_few_phases(): + with pytest.raises(ValueError, match="at least 3 phases"): + to_dq(np.zeros((len(THETA), 2)), THETA, 1) From 2ddad108f53ab24516ec98b609f2b06a20ba40f4 Mon Sep 17 00:00:00 2001 From: Jan Laksar Date: Tue, 15 Sep 2026 08:07:19 +0200 Subject: [PATCH 2/2] Independent creation of Output variables separate files applicable to different number of phases, excitation (current vs. voltage), rotation direction, names of phases, ... --- Ansys_transforms/mech_offset.py | 39 +++ Ansys_transforms/output_vars.py | 100 +++++++ Ansys_transforms/run.py | 27 ++ Ansys_transforms/symbolic_builder.py | 382 +++++++++++++++++++++++++++ Ansys_transforms/user_config.py | 114 ++++++++ Ansys_transforms/validate_inputs.py | 278 +++++++++++++++++++ Ansys_transforms/write.py | 68 +++++ 7 files changed, 1008 insertions(+) create mode 100644 Ansys_transforms/mech_offset.py create mode 100644 Ansys_transforms/output_vars.py create mode 100644 Ansys_transforms/run.py create mode 100644 Ansys_transforms/symbolic_builder.py create mode 100644 Ansys_transforms/user_config.py create mode 100644 Ansys_transforms/validate_inputs.py create mode 100644 Ansys_transforms/write.py diff --git a/Ansys_transforms/mech_offset.py b/Ansys_transforms/mech_offset.py new file mode 100644 index 0000000..a37d675 --- /dev/null +++ b/Ansys_transforms/mech_offset.py @@ -0,0 +1,39 @@ +import math + + +def compute_mech_offset_from_winding(inputs): + + Q = inputs.params["slot_number"]# number of slots + p = inputs.params["pole_pairs"] # pole pairs + m = inputs.params["num_phases"] # phases + u2 = inputs.params["layers"] # 1 or 2 + ys = inputs.params["coil_pitch"] # coil pitch (slots) + + alpha_s = p * 360.0 / Q # electrical angle between slots + phase_sector = 180.0 / m + + angles = [] + + # slot iteration step. for a single-layer winding, the front coil sides are only in odd slots + step = 1 if u2 == 2 else 2 + + for k in range(1, Q + 1, step): + + theta = ((k - 1) * alpha_s) % 180 + + # slot belongs to phase A sector + if 0 <= theta < phase_sector: + angles.append(theta) + + if not angles: + raise ValueError("No slots found for phase A. Check winding parameters.") + + theta_avg = sum(angles) / len(angles) + + # shift from slot side to coil axis + theta_shift = (ys / 2.0) * alpha_s + + phAaxis = theta_avg + theta_shift + alpha_s/2 #in FEM model, the angle 0 deg is between first and last slots + mech_offset = inputs.params["rotation_sign"]*(phAaxis - 90.0)/p # the north PM location is at 90 el. degrees. For the negative rotation, the stator winding is built inversely by Motor-CAD + inputs.params["mechanical_offset_deg"] = mech_offset + return inputs \ No newline at end of file diff --git a/Ansys_transforms/output_vars.py b/Ansys_transforms/output_vars.py new file mode 100644 index 0000000..8dd44a6 --- /dev/null +++ b/Ansys_transforms/output_vars.py @@ -0,0 +1,100 @@ +output_vars = { + 'PolePairs': "2", + 'RotSign': "1", + 'MechOffsetDeg': "0.0", + 'Rs': "19", + 'Lew': "0", + 'theta_el': "RotSign*(Moving1.Position - MechOffsetDeg/180*pi) * PolePairs - pi", + 'cos0_1': "cos(1*(theta_el - 2*PI*0/5))", + 'sin0_1': "sin(-1*(theta_el - 2*PI*0/5))", + 'cos1_1': "cos(1*(theta_el - 2*PI*1/5))", + 'sin1_1': "sin(-1*(theta_el - 2*PI*1/5))", + 'cos2_1': "cos(1*(theta_el - 2*PI*2/5))", + 'sin2_1': "sin(-1*(theta_el - 2*PI*2/5))", + 'cos3_1': "cos(1*(theta_el - 2*PI*3/5))", + 'sin3_1': "sin(-1*(theta_el - 2*PI*3/5))", + 'cos4_1': "cos(1*(theta_el - 2*PI*4/5))", + 'sin4_1': "sin(-1*(theta_el - 2*PI*4/5))", + 'cos0_3': "cos(3*(theta_el - 2*PI*0/5))", + 'sin0_3': "sin(-3*(theta_el - 2*PI*0/5))", + 'cos1_3': "cos(3*(theta_el - 2*PI*1/5))", + 'sin1_3': "sin(-3*(theta_el - 2*PI*1/5))", + 'cos2_3': "cos(3*(theta_el - 2*PI*2/5))", + 'sin2_3': "sin(-3*(theta_el - 2*PI*2/5))", + 'cos3_3': "cos(3*(theta_el - 2*PI*3/5))", + 'sin3_3': "sin(-3*(theta_el - 2*PI*3/5))", + 'cos4_3': "cos(3*(theta_el - 2*PI*4/5))", + 'sin4_3': "sin(-3*(theta_el - 2*PI*4/5))", + 'Flux_d1': "(FluxLinkage(PhaseA)*cos0_1 + FluxLinkage(PhaseB)*cos1_1 + FluxLinkage(PhaseC)*cos2_1 + FluxLinkage(PhaseD)*cos3_1 + FluxLinkage(PhaseE)*cos4_1) * 2/5", + 'Flux_q1': "(FluxLinkage(PhaseA)*sin0_1 + FluxLinkage(PhaseB)*sin1_1 + FluxLinkage(PhaseC)*sin2_1 + FluxLinkage(PhaseD)*sin3_1 + FluxLinkage(PhaseE)*sin4_1) * 2/5", + 'Flux_d3': "(FluxLinkage(PhaseA)*cos0_3 + FluxLinkage(PhaseB)*cos1_3 + FluxLinkage(PhaseC)*cos2_3 + FluxLinkage(PhaseD)*cos3_3 + FluxLinkage(PhaseE)*cos4_3) * 2/5", + 'Flux_q3': "(FluxLinkage(PhaseA)*sin0_3 + FluxLinkage(PhaseB)*sin1_3 + FluxLinkage(PhaseC)*sin2_3 + FluxLinkage(PhaseD)*sin3_3 + FluxLinkage(PhaseE)*sin4_3) * 2/5", + 'Vind_d1': "(InducedVoltage(PhaseA)*cos0_1 + InducedVoltage(PhaseB)*cos1_1 + InducedVoltage(PhaseC)*cos2_1 + InducedVoltage(PhaseD)*cos3_1 + InducedVoltage(PhaseE)*cos4_1) * 2/5", + 'Vind_q1': "(InducedVoltage(PhaseA)*sin0_1 + InducedVoltage(PhaseB)*sin1_1 + InducedVoltage(PhaseC)*sin2_1 + InducedVoltage(PhaseD)*sin3_1 + InducedVoltage(PhaseE)*sin4_1) * 2/5", + 'Vind_d3': "(InducedVoltage(PhaseA)*cos0_3 + InducedVoltage(PhaseB)*cos1_3 + InducedVoltage(PhaseC)*cos2_3 + InducedVoltage(PhaseD)*cos3_3 + InducedVoltage(PhaseE)*cos4_3) * 2/5", + 'Vind_q3': "(InducedVoltage(PhaseA)*sin0_3 + InducedVoltage(PhaseB)*sin1_3 + InducedVoltage(PhaseC)*sin2_3 + InducedVoltage(PhaseD)*sin3_3 + InducedVoltage(PhaseE)*sin4_3) * 2/5", + 'V_A': "InducedVoltage(PhaseA) + Rs*InputCurrent(PhaseA) + Lew*ddt(InputCurrent(PhaseA))", + 'V_B': "InducedVoltage(PhaseB) + Rs*InputCurrent(PhaseB) + Lew*ddt(InputCurrent(PhaseB))", + 'V_C': "InducedVoltage(PhaseC) + Rs*InputCurrent(PhaseC) + Lew*ddt(InputCurrent(PhaseC))", + 'V_D': "InducedVoltage(PhaseD) + Rs*InputCurrent(PhaseD) + Lew*ddt(InputCurrent(PhaseD))", + 'V_E': "InducedVoltage(PhaseE) + Rs*InputCurrent(PhaseE) + Lew*ddt(InputCurrent(PhaseE))", + 'V_AC': "V_A - V_C", + 'V_BD': "V_B - V_D", + 'V_CE': "V_C - V_E", + 'V_DA': "V_D - V_A", + 'V_EB': "V_E - V_B", + 'Vterm_A': "1/5*(2*V_AC + -1*V_BD + 1*V_CE + -2*V_DA)", + 'Vterm_B': "1/5*(2*V_AC + 4*V_BD + 1*V_CE + 3*V_DA)", + 'Vterm_C': "1/5*(-3*V_AC + -1*V_BD + 1*V_CE + -2*V_DA)", + 'Vterm_D': "1/5*(2*V_AC + -1*V_BD + 1*V_CE + 3*V_DA)", + 'Vterm_E': "1/5*(-3*V_AC + -1*V_BD + -4*V_CE + -2*V_DA)", + 'I_d1': "(InputCurrent(PhaseA)*cos0_1 + InputCurrent(PhaseB)*cos1_1 + InputCurrent(PhaseC)*cos2_1 + InputCurrent(PhaseD)*cos3_1 + InputCurrent(PhaseE)*cos4_1) * 2/5", + 'I_q1': "(InputCurrent(PhaseA)*sin0_1 + InputCurrent(PhaseB)*sin1_1 + InputCurrent(PhaseC)*sin2_1 + InputCurrent(PhaseD)*sin3_1 + InputCurrent(PhaseE)*sin4_1) * 2/5", + 'I_d3': "(InputCurrent(PhaseA)*cos0_3 + InputCurrent(PhaseB)*cos1_3 + InputCurrent(PhaseC)*cos2_3 + InputCurrent(PhaseD)*cos3_3 + InputCurrent(PhaseE)*cos4_3) * 2/5", + 'I_q3': "(InputCurrent(PhaseA)*sin0_3 + InputCurrent(PhaseB)*sin1_3 + InputCurrent(PhaseC)*sin2_3 + InputCurrent(PhaseD)*sin3_3 + InputCurrent(PhaseE)*sin4_3) * 2/5", + 'V_d1': "(V_A*cos0_1 + V_B*cos1_1 + V_C*cos2_1 + V_D*cos3_1 + V_E*cos4_1) * 2/5", + 'V_q1': "(V_A*sin0_1 + V_B*sin1_1 + V_C*sin2_1 + V_D*sin3_1 + V_E*sin4_1) * 2/5", + 'V_d3': "(V_A*cos0_3 + V_B*cos1_3 + V_C*cos2_3 + V_D*cos3_3 + V_E*cos4_3) * 2/5", + 'V_q3': "(V_A*sin0_3 + V_B*sin1_3 + V_C*sin2_3 + V_D*sin3_3 + V_E*sin4_3) * 2/5", + 'L0d_1': "L(PhaseA,PhaseA)*cos0_1 + L(PhaseA,PhaseB)*cos1_1 + L(PhaseA,PhaseC)*cos2_1 + L(PhaseA,PhaseD)*cos3_1 + L(PhaseA,PhaseE)*cos4_1", + 'L0q_1': "L(PhaseA,PhaseA)*sin0_1 + L(PhaseA,PhaseB)*sin1_1 + L(PhaseA,PhaseC)*sin2_1 + L(PhaseA,PhaseD)*sin3_1 + L(PhaseA,PhaseE)*sin4_1", + 'L1d_1': "L(PhaseB,PhaseA)*cos0_1 + L(PhaseB,PhaseB)*cos1_1 + L(PhaseB,PhaseC)*cos2_1 + L(PhaseB,PhaseD)*cos3_1 + L(PhaseB,PhaseE)*cos4_1", + 'L1q_1': "L(PhaseB,PhaseA)*sin0_1 + L(PhaseB,PhaseB)*sin1_1 + L(PhaseB,PhaseC)*sin2_1 + L(PhaseB,PhaseD)*sin3_1 + L(PhaseB,PhaseE)*sin4_1", + 'L2d_1': "L(PhaseC,PhaseA)*cos0_1 + L(PhaseC,PhaseB)*cos1_1 + L(PhaseC,PhaseC)*cos2_1 + L(PhaseC,PhaseD)*cos3_1 + L(PhaseC,PhaseE)*cos4_1", + 'L2q_1': "L(PhaseC,PhaseA)*sin0_1 + L(PhaseC,PhaseB)*sin1_1 + L(PhaseC,PhaseC)*sin2_1 + L(PhaseC,PhaseD)*sin3_1 + L(PhaseC,PhaseE)*sin4_1", + 'L3d_1': "L(PhaseD,PhaseA)*cos0_1 + L(PhaseD,PhaseB)*cos1_1 + L(PhaseD,PhaseC)*cos2_1 + L(PhaseD,PhaseD)*cos3_1 + L(PhaseD,PhaseE)*cos4_1", + 'L3q_1': "L(PhaseD,PhaseA)*sin0_1 + L(PhaseD,PhaseB)*sin1_1 + L(PhaseD,PhaseC)*sin2_1 + L(PhaseD,PhaseD)*sin3_1 + L(PhaseD,PhaseE)*sin4_1", + 'L4d_1': "L(PhaseE,PhaseA)*cos0_1 + L(PhaseE,PhaseB)*cos1_1 + L(PhaseE,PhaseC)*cos2_1 + L(PhaseE,PhaseD)*cos3_1 + L(PhaseE,PhaseE)*cos4_1", + 'L4q_1': "L(PhaseE,PhaseA)*sin0_1 + L(PhaseE,PhaseB)*sin1_1 + L(PhaseE,PhaseC)*sin2_1 + L(PhaseE,PhaseD)*sin3_1 + L(PhaseE,PhaseE)*sin4_1", + 'L0d_3': "L(PhaseA,PhaseA)*cos0_3 + L(PhaseA,PhaseB)*cos1_3 + L(PhaseA,PhaseC)*cos2_3 + L(PhaseA,PhaseD)*cos3_3 + L(PhaseA,PhaseE)*cos4_3", + 'L0q_3': "L(PhaseA,PhaseA)*sin0_3 + L(PhaseA,PhaseB)*sin1_3 + L(PhaseA,PhaseC)*sin2_3 + L(PhaseA,PhaseD)*sin3_3 + L(PhaseA,PhaseE)*sin4_3", + 'L1d_3': "L(PhaseB,PhaseA)*cos0_3 + L(PhaseB,PhaseB)*cos1_3 + L(PhaseB,PhaseC)*cos2_3 + L(PhaseB,PhaseD)*cos3_3 + L(PhaseB,PhaseE)*cos4_3", + 'L1q_3': "L(PhaseB,PhaseA)*sin0_3 + L(PhaseB,PhaseB)*sin1_3 + L(PhaseB,PhaseC)*sin2_3 + L(PhaseB,PhaseD)*sin3_3 + L(PhaseB,PhaseE)*sin4_3", + 'L2d_3': "L(PhaseC,PhaseA)*cos0_3 + L(PhaseC,PhaseB)*cos1_3 + L(PhaseC,PhaseC)*cos2_3 + L(PhaseC,PhaseD)*cos3_3 + L(PhaseC,PhaseE)*cos4_3", + 'L2q_3': "L(PhaseC,PhaseA)*sin0_3 + L(PhaseC,PhaseB)*sin1_3 + L(PhaseC,PhaseC)*sin2_3 + L(PhaseC,PhaseD)*sin3_3 + L(PhaseC,PhaseE)*sin4_3", + 'L3d_3': "L(PhaseD,PhaseA)*cos0_3 + L(PhaseD,PhaseB)*cos1_3 + L(PhaseD,PhaseC)*cos2_3 + L(PhaseD,PhaseD)*cos3_3 + L(PhaseD,PhaseE)*cos4_3", + 'L3q_3': "L(PhaseD,PhaseA)*sin0_3 + L(PhaseD,PhaseB)*sin1_3 + L(PhaseD,PhaseC)*sin2_3 + L(PhaseD,PhaseD)*sin3_3 + L(PhaseD,PhaseE)*sin4_3", + 'L4d_3': "L(PhaseE,PhaseA)*cos0_3 + L(PhaseE,PhaseB)*cos1_3 + L(PhaseE,PhaseC)*cos2_3 + L(PhaseE,PhaseD)*cos3_3 + L(PhaseE,PhaseE)*cos4_3", + 'L4q_3': "L(PhaseE,PhaseA)*sin0_3 + L(PhaseE,PhaseB)*sin1_3 + L(PhaseE,PhaseC)*sin2_3 + L(PhaseE,PhaseD)*sin3_3 + L(PhaseE,PhaseE)*sin4_3", + 'Ld1': "(L0d_1*cos0_1 + L1d_1*cos1_1 + L2d_1*cos2_1 + L3d_1*cos3_1 + L4d_1*cos4_1) * 2/5", + 'Ld1q1': "(L0d_1*sin0_1 + L1d_1*sin1_1 + L2d_1*sin2_1 + L3d_1*sin3_1 + L4d_1*sin4_1) * 2/5", + 'Lq1d1': "(L0q_1*cos0_1 + L1q_1*cos1_1 + L2q_1*cos2_1 + L3q_1*cos3_1 + L4q_1*cos4_1) * 2/5", + 'Lq1': "(L0q_1*sin0_1 + L1q_1*sin1_1 + L2q_1*sin2_1 + L3q_1*sin3_1 + L4q_1*sin4_1) * 2/5", + 'Ld1d3': "(L0d_1*cos0_3 + L1d_1*cos1_3 + L2d_1*cos2_3 + L3d_1*cos3_3 + L4d_1*cos4_3) * 2/5", + 'Ld1q3': "(L0d_1*sin0_3 + L1d_1*sin1_3 + L2d_1*sin2_3 + L3d_1*sin3_3 + L4d_1*sin4_3) * 2/5", + 'Lq1d3': "(L0q_1*cos0_3 + L1q_1*cos1_3 + L2q_1*cos2_3 + L3q_1*cos3_3 + L4q_1*cos4_3) * 2/5", + 'Lq1q3': "(L0q_1*sin0_3 + L1q_1*sin1_3 + L2q_1*sin2_3 + L3q_1*sin3_3 + L4q_1*sin4_3) * 2/5", + 'Ld3d1': "(L0d_3*cos0_1 + L1d_3*cos1_1 + L2d_3*cos2_1 + L3d_3*cos3_1 + L4d_3*cos4_1) * 2/5", + 'Ld3q1': "(L0d_3*sin0_1 + L1d_3*sin1_1 + L2d_3*sin2_1 + L3d_3*sin3_1 + L4d_3*sin4_1) * 2/5", + 'Lq3d1': "(L0q_3*cos0_1 + L1q_3*cos1_1 + L2q_3*cos2_1 + L3q_3*cos3_1 + L4q_3*cos4_1) * 2/5", + 'Lq3q1': "(L0q_3*sin0_1 + L1q_3*sin1_1 + L2q_3*sin2_1 + L3q_3*sin3_1 + L4q_3*sin4_1) * 2/5", + 'Ld3': "(L0d_3*cos0_3 + L1d_3*cos1_3 + L2d_3*cos2_3 + L3d_3*cos3_3 + L4d_3*cos4_3) * 2/5", + 'Ld3q3': "(L0d_3*sin0_3 + L1d_3*sin1_3 + L2d_3*sin2_3 + L3d_3*sin3_3 + L4d_3*sin4_3) * 2/5", + 'Lq3d3': "(L0q_3*cos0_3 + L1q_3*cos1_3 + L2q_3*cos2_3 + L3q_3*cos3_3 + L4q_3*cos4_3) * 2/5", + 'Lq3': "(L0q_3*sin0_3 + L1q_3*sin1_3 + L2q_3*sin2_3 + L3q_3*sin3_3 + L4q_3*sin4_3) * 2/5", + 'Flux_e_d1': "Flux_d1 - (Ld1*I_d1 + Ld1q1*I_q1 + Ld1d3*I_d3 + Ld1q3*I_q3)", + 'Flux_e_q1': "Flux_q1 - (Lq1d1*I_d1 + Lq1*I_q1 + Lq1d3*I_d3 + Lq1q3*I_q3)", + 'Flux_e_d3': "Flux_d3 - (Ld3d1*I_d1 + Ld3q1*I_q1 + Ld3*I_d3 + Ld3q3*I_q3)", + 'Flux_e_q3': "Flux_q3 - (Lq3d1*I_d1 + Lq3q1*I_q1 + Lq3d3*I_d3 + Lq3*I_q3)", + 'Torque_dq': "5/2*PolePairs*(1*(Flux_d1*I_q1 - Flux_q1*I_d1) + 3*(Flux_d3*I_q3 - Flux_q3*I_d3))", +} diff --git a/Ansys_transforms/run.py b/Ansys_transforms/run.py new file mode 100644 index 0000000..9974664 --- /dev/null +++ b/Ansys_transforms/run.py @@ -0,0 +1,27 @@ +from validate_inputs import AoutvarConfig +from symbolic_builder import SymbolicBuilder +from mech_offset import compute_mech_offset_from_winding +from write import write_aoutvar + +def run(CONFIG): + + print("\n======================================") + print(" AOUTVAR GENERATOR STARTED") + print("======================================") + + # 1) Validate configuration + inputs = AoutvarConfig(CONFIG) + + # 2) Automatically compute the electrical offset of the winding - if chosen + if inputs.params["mechanical_offset_deg"] == "auto": + inputs = compute_mech_offset_from_winding(inputs) + + # 2) Build symbolic variables + builder = SymbolicBuilder(inputs.params) + + # 3) Write output file + write_aoutvar( builder.variables, inputs.params) + + print("\n[OK] output.aoutvar successfully created.") + print("======================================\n") + diff --git a/Ansys_transforms/symbolic_builder.py b/Ansys_transforms/symbolic_builder.py new file mode 100644 index 0000000..5ab71e0 --- /dev/null +++ b/Ansys_transforms/symbolic_builder.py @@ -0,0 +1,382 @@ +""" +SYMBOLIC dq EXPRESSION GENERATOR +============================================================ + +Generates symbolic dq transformation expressions for +multiphase machines. + +The expressions are written in a format compatible with +ANSYS Electronics Desktop (.aoutvar). + +Supported features +------------------ + +• arbitrary number of phases +• harmonic dq subspaces +• dq projection of currents, voltages and flux linkages +• dq inductance matrices +• excitation flux calculation + +All generated expressions are stored as strings and later +written to the .aoutvar file. + +------------------------------------------------------------ +Author: Jan Laksar +------------------------------------------------------------ +""" + +# ============================================================ +# Imports +# ============================================================ + +# Third-party libraries +import numpy as np + +# ============================================================ +# Symbolic builder class +# ============================================================ + +class SymbolicBuilder: + + # ------------------------------------------------------------ + # Initialization + # ------------------------------------------------------------ + + def __init__(self, inputs): + self.inputs = inputs + self.variables = {} + + self._build_transform_constants() + self._build_theta() + self._build_trigonometry() + self._build_dq_projection("Flux", "FluxLinkage", "Wb") + self._build_dq_projection("Vind", "InducedVoltage", "V") + self._build_voltages() + if self.inputs["excitation_type"] == "voltage": + self._build_dq_projection("I", "Current", "A") + self._build_dq_projection("V", "InputVoltage", "V") + else: + self._build_dq_projection("I", "InputCurrent", "A") + self._build_dq_voltages() + if self.inputs["include_inductances"]==True: + self._build_phase_inductance_projection() + self._build_dq_inductances() + self._build_excitation_fluxes() + self._build_electromagnetic_torque() + + print("\n--- Checking configuration ---\n") + + # ------------------------------------------------------------ + # Transformation constants + # ------------------------------------------------------------ + + def _build_transform_constants(self): + cfg = self.inputs + self.add("PolePairs",str(cfg["pole_pairs"])) + self.add("RotSign",str(cfg["rotation_sign"])) + self.add("MechOffsetDeg",str(cfg["mechanical_offset_deg"])) + self.add("Rs",str(cfg["resistance_stat"])) + self.add("Lew",str(cfg["inductance_endw"])) + + # ------------------------------------------------------------ + # Electrical rotor angle + # ------------------------------------------------------------ + + def _build_theta(self): + motion = self.inputs["motion_setup"] + expr = f"RotSign*({motion}.Position - MechOffsetDeg/180*pi) * PolePairs - pi" + self.add("theta_el",expr,False,"deg","Uses20Log") + + # ------------------------------------------------------------ + # Trigonometric basis functions + # ------------------------------------------------------------ + + def _build_trigonometry(self): + m = self.inputs["num_phases"] + harmonics = self.inputs["harmonic_orders"] + + for mu in harmonics: + for k in range(m): + + name_cos = f"cos{k}_{mu}" + name_sin = f"sin{k}_{mu}" + angle = f"{mu}*(theta_el - 2*PI*{k}/{m})" + + self.add(name_cos, f"cos({angle})", False) + self.add(name_sin, f"sin(-{angle})", False) + + # ------------------------------------------------------------------- + # dq projection of phase currents, flux linkages and induced voltages + # ------------------------------------------------------------------- + + def _build_dq_projection(self, name, ansys_func, unit): + m = self.inputs["num_phases"] + phases = self.inputs["phase_names"] + harmonics = self.inputs["harmonic_orders"] + + for mu in harmonics: + terms_d = [] + terms_q = [] + for k, phase in enumerate(phases): + name_cos = f"cos{k}_{mu}" + name_sin = f"sin{k}_{mu}" + terms_d.append(f"{ansys_func}({phase})*{name_cos}") + terms_q.append(f"{ansys_func}({phase})*{name_sin}") + + expr_d = "(" + " + ".join(terms_d) + f") * 2/{m}" + expr_q = "(" + " + ".join(terms_q) + f") * 2/{m}" + + self.add(f"{name}_d{mu}", expr_d, False, unit, "Uses20Log") + self.add(f"{name}_q{mu}", expr_q, False, unit, "Uses20Log") + + # ------------------------------------------------------------ + # Projection of phase inductances + # ------------------------------------------------------------ + + def _build_phase_inductance_projection(self): + phases = self.inputs["phase_names"] + harmonics = self.inputs["harmonic_orders"] + + for mu in harmonics: + for i, phase_i in enumerate(phases): + terms_d = [] + terms_q = [] + name_d = f"L{i}d_{mu}" + name_q = f"L{i}q_{mu}" + + for k, phase_k in enumerate(phases): + name_cos = f"cos{k}_{mu}" + name_sin = f"sin{k}_{mu}" + + terms_d.append(f"L({phase_i},{phase_k})*{name_cos}") + terms_q.append(f"L({phase_i},{phase_k})*{name_sin}") + + expr_d = f" + ".join(terms_d) + expr_q = f" + ".join(terms_q) + + self.add(name_d, expr_d, False, "nH", "Uses20Log") + self.add(name_q, expr_q, False, "nH", "Uses20Log") + + # ------------------------------------------------------------ + # dq inductance matrices + # ------------------------------------------------------------ + + def _build_dq_inductances(self): + m = self.inputs["num_phases"] + phases = self.inputs["phase_names"] + harmonics = self.inputs["harmonic_orders"] + + for mu in harmonics: + for nu in harmonics: + terms_dd = [] + terms_dq = [] + terms_qq = [] + terms_qd = [] + + if mu==nu: + name_dd = f"Ld{mu}" + name_qq = f"Lq{mu}" + else: + name_dd = f"Ld{mu}d{nu}" + name_qq = f"Lq{mu}q{nu}" + name_dq = f"Ld{mu}q{nu}" + name_qd = f"Lq{mu}d{nu}" + + for k in range(len(phases)): + + coef_d = f"L{k}d_{mu}" + coef_q = f"L{k}q_{mu}" + + name_cos = f"cos{k}_{nu}" + name_sin = f"sin{k}_{nu}" + + # Ldd = dΨdμ / didν + terms_dd.append(f"{coef_d}*{name_cos}") + # Ldq = dΨdμ / diqν + terms_dq.append(f"{coef_d}*{name_sin}") + # Lqq = dΨqμ / diqν + terms_qq.append(f"{coef_q}*{name_sin}") + # Lqd = dΨqμ / didν + terms_qd.append(f"{coef_q}*{name_cos}") + + if self.inputs["inductance_mode"] == "full": + expr_dd = f"(" + " + ".join(terms_dd) + f") * 2/{m}" + expr_dq = f"(" + " + ".join(terms_dq) + f") * 2/{m}" + expr_qd = f"(" + " + ".join(terms_qd) + f") * 2/{m}" + expr_qq = f"(" + " + ".join(terms_qq) + f") * 2/{m}" + self.add(name_dd, expr_dd, False, "nH", "Uses20Log") + self.add(name_dq, expr_dq, False, "nH", "Uses20Log") + self.add(name_qd, expr_qd, False, "nH", "Uses20Log") + self.add(name_qq, expr_qq, False, "nH", "Uses20Log") + elif mu==nu: + expr_dd = f"(" + " + ".join(terms_dd) + f") * 2/{m}" + expr_qq = f"(" + " + ".join(terms_qq) + f") * 2/{m}" + self.add(name_dd, expr_dd, False, "nH", "Uses20Log") + self.add(name_qq, expr_qq, False, "nH", "Uses20Log") + + # ------------------------------------------------------------ + # Excitation flux calculation + # ------------------------------------------------------------ + + def _build_excitation_fluxes(self): + harmonics = self.inputs["harmonic_orders"] + + for mu in harmonics: + terms_d = [] + terms_q = [] + + for nu in harmonics: + if mu==nu: + name_dd = f"Ld{mu}" + name_qq = f"Lq{mu}" + else: + name_dd = f"Ld{mu}d{nu}" + name_qq = f"Lq{mu}q{nu}" + name_dq = f"Ld{mu}q{nu}" + name_qd = f"Lq{mu}d{nu}" + + if self.inputs["inductance_mode"]=="full": + terms_d.append(f"{name_dd}*I_d{nu}") + terms_d.append(f"{name_dq}*I_q{nu}") + terms_q.append(f"{name_qd}*I_d{nu}") + terms_q.append(f"{name_qq}*I_q{nu}") + elif mu==nu: + terms_d.append(f"{name_dd}*I_d{nu}") + terms_q.append(f"{name_qq}*I_q{nu}") + + expr_d = " + ".join(terms_d) + expr_q = " + ".join(terms_q) + + self.add(f"Flux_e_d{mu}",f"Flux_d{mu} - ({expr_d})", False, "''", "Uses20Log") + self.add(f"Flux_e_q{mu}",f"Flux_q{mu} - ({expr_q})", False, "''", "Uses20Log") + + # ------------------------------------------------------------ + # Phase and terminal voltages + # ------------------------------------------------------------ + + def _build_voltages(self): + m = self.inputs["num_phases"] + phases = self.inputs["phase_names"] + codes = self.inputs["phase_codes"] + Dm = _build_phase_from_line_matrix(m) + + if m % 2 == 0: + delta = m // 2 + else: + delta = (m - 1) // 2 + + phase_vars = [] + line_vars = [] + term_vars = [] + + for k, phase_a in enumerate(phases): + phase_b = phases[(k + delta) % m] + code_a = codes[k] + code_b = codes[(k + delta) % m] + + var_name_ll = f"V_{code_a}{code_b}" + if self.inputs["excitation_type"] == "voltage": + expr_ll = f"InputVoltage({phase_a}) - InputVoltage({phase_b})" + else: + var_name_ph = f"V_{code_a}" + var_name_term = f"Vterm_{code_a}" + expr_ph = f"InducedVoltage({phase_a}) + Rs*InputCurrent({phase_a}) + Lew*deriv(InputCurrent({phase_a}))" + expr_ll = f"V_{code_a} - V_{code_b}" + + terms = [] + for i in range(m-1): + ph_a = codes[i] + ph_b = codes[(i + delta) % m] + + var_ll = f"V_{ph_a}{ph_b}" + + terms.append(f"{Dm[k,i]}*{var_ll}") + + expr_term = " + ".join(terms) + expr_term = f"1/{m}*({expr_term})" + + phase_vars.append((var_name_ph, expr_ph)) + term_vars.append((var_name_term, expr_term)) + + line_vars.append((var_name_ll, expr_ll)) + + for name, expr in phase_vars: + self.add(name, expr, False, "''", "Uses20Log") + for name, expr in line_vars: + self.add(name, expr, False, "''", "Uses20Log") + for name, expr in term_vars: + self.add(name, expr, False, "''", "Uses20Log") + + # ------------------------------------------------------------ + # dq projection of phase voltages for current excitation + # ------------------------------------------------------------ + + def _build_dq_voltages(self): + m = self.inputs["num_phases"] + phases = self.inputs["phase_names"] + codes = self.inputs["phase_codes"] + harmonics = self.inputs["harmonic_orders"] + + for mu in harmonics: + terms_d = [] + terms_q = [] + for k, code in enumerate(codes): + name_cos = f"cos{k}_{mu}" + name_sin = f"sin{k}_{mu}" + terms_d.append(f"V_{code}*{name_cos}") + terms_q.append(f"V_{code}*{name_sin}") + + expr_d = "(" + " + ".join(terms_d) + f") * 2/{m}" + expr_q = "(" + " + ".join(terms_q) + f") * 2/{m}" + + self.add(f"V_d{mu}", expr_d, False, "''", "Uses20Log") + self.add(f"V_q{mu}", expr_q, False, "''", "Uses20Log") + + # ------------------------------------------------------------ + # dq projection of electromagnetic torque + # ------------------------------------------------------------ + + def _build_electromagnetic_torque(self): + m = self.inputs["num_phases"] + harmonics = self.inputs["harmonic_orders"] + + terms = [] + + for mu in harmonics: + term = f"{mu}*(Flux_d{mu}*I_q{mu} - Flux_q{mu}*I_d{mu})" + terms.append(term) + + expr = f"{m}/2*PolePairs*(" + " + ".join(terms) + ")" + + self.add(f"Torque_dq", expr, False, "''", "Uses20Log") + + # ------------------------------------------------------------ + # Variable registration helper + # ------------------------------------------------------------ + + def add(self, name, expression, is_constant=True, unit="''", db_flag="dBTypeDoesntCare"): + self.variables[name] = { + "expression": expression, + "unit": unit, + "db_flag": db_flag, + "is_constant": is_constant + } + +# ============================================================ +# Phase-to-line voltage transformation matrix +# ============================================================ + +def _build_phase_from_line_matrix(m): + if m % 2 == 0: + k = m//2 + else: + k = (m-1)//2 + + D = np.zeros((m,m)) + + for i in range(m): + D[i,i] = 1 + D[i,(i+k)%m] = -1 + D[-1,:] = 1 + Dinv = np.linalg.inv(D) + return np.round(Dinv * m).astype(int) \ No newline at end of file diff --git a/Ansys_transforms/user_config.py b/Ansys_transforms/user_config.py new file mode 100644 index 0000000..70e768b --- /dev/null +++ b/Ansys_transforms/user_config.py @@ -0,0 +1,114 @@ +""" +AOUTVAR GENERATOR +----------------- +USER CONFIGURATION FILE +----------------------- +Edit CONFIG below and run this file. + +All unspecified parameters will use safe defaults. +See README for detailed explanation of each parameter. +""" + +CONFIG = { + # 1) MACHINE DEFINITION + + # Number of pole pairs (mandatory): Integer + "pole_pairs": 3, + # Number of phases (mandatory): Integer + "num_phases": 5, + + # 2) MECHANICAL SETTINGS + + # Name of motion setup in Ansys + # Default: "Moving1" + "motion_setup": "Moving1", + + # Rotation direction: + # +1 = positive (CCW) + # -1 = negative (CW) + # Default: +1 + "rotation_sign": +1, + + + # Additional mechanical offset in degrees + # "auto" - automatically calculates when number of slots, winding layers and winding step are defined + # Default: 0.0 + "mechanical_offset_deg": 0.0, + + # Additional winding information for the automatic mechanical offset calcualtion + # layers == number of coil sides per slot: 1 or 2 + # coil pitch in number of slots: MUST be odd number of a single-layer winding + # Default: 0 + "slot_number": 12, + "layers": 2, + "coil_pitch": 1, + + # 3) EXCITATION TYPE + # "Current" (default) + # "Voltage" + "excitation_type": "Voltage", + + # 4) PHASE NAMING + + # Base name for automatic phase naming (case sensitive) + # Default: "Phase" + "phase_base_name": "Phase", + + # Index style: + # "upper" → PhaseA, PhaseB, PhaseC + # "lower" → Phasea, Phaseb, Phasec + # "number" → Phase1, Phase2, Phase3 + # Default: "upper" + "phase_index_style": "upper", + + # Optional: explicit phase names + # Example: + # ["U1","V1","W1","U2","V2","W2"] + # None + # If provided, this overrides base name and index style. + # "custom_phase_names": ["WG_Ph1_P1","WG_Ph4_P1","WG_Ph7_P1","WG_Ph2_P1","WG_Ph5_P1","WG_Ph8_P1","WG_Ph3_P1","WG_Ph6_P1","WG_Ph9_P1"], + "custom_phase_names": None, + + # 5) HARMONIC SUBSPACES + # Spatial harmonic orders ν + # Example: + # [1] → fundamental only + # [1, 3] → fundamental + 3rd harmonic subspace + # [1, 5, 7] → extended harmonic analysis + # Default: [1] + "harmonic_orders": [1, 3], + + # 6) INDUCTANCE PROJECTION + + # Generate inductance expressions (Ld, Lq) + # Default: True + # False + "include_inductances": True, + + # Inductance mode: + # "full" → full projection including cross-coupling + # "main" → only main diagonal Ld and Lq + # Default: "full" + "inductance_mode": "main", + + # 7) OPTIONALS - values which can be used in the "Current" Excitation type to calculate the terminal voltages + # Stator winding resistance [Ohm] + # Default: 0 + "resistance_stat": 0.0, + # Endwinding inductance [H] + # Default: 0 + "inductance_endw": 0.0e-6, +} + + +# ============================================================ +# DO NOT MODIFY BELOW THIS LINE +# ============================================================ + +from run import run + +def main(): + run(CONFIG) + +if __name__ == "__main__": + main() diff --git a/Ansys_transforms/validate_inputs.py b/Ansys_transforms/validate_inputs.py new file mode 100644 index 0000000..4175ee7 --- /dev/null +++ b/Ansys_transforms/validate_inputs.py @@ -0,0 +1,278 @@ +""" +CONFIGURATION VALIDATION MODULE +============================================================ + +Responsible for: + +• validating user configuration parameters +• filling missing parameters with default values +• normalizing input formats +• preparing derived parameters used by the generator + +This module performs NO symbolic generation. + +Output: + validated configuration stored in `self.params` + +------------------------------------------------------------ +Author: Jan Laksar +------------------------------------------------------------ +""" + +# ============================================================ +# Configuration object +# ============================================================ + +class AoutvarConfig: + + # ============================================================ + # Default parameter values + # ============================================================ + + DEFAULTS = { + # Machine + "pole_pairs": 2, + "num_phases": 3, + # Mechanical + "motion_setup": "Moving1", + "rotation_sign": 1, + "mechanical_offset_deg": 0.0, + "slot_number": 0, + "layers": 0, + "coil_pitch": 0, + # Excitation + "excitation_type": "current", + # Phase system + "phase_base_name": "Phase", + "phase_index_style": "upper", + "custom_phase_names": None, + + # Harmonics + "harmonic_orders": [1], + + # Inductances + "include_inductances": True, + "inductance_mode": "full", + + # Winding parameters + } + + # ============================================================ + # Validation rules + # ============================================================ + + VALIDATION_RULES = { + "pole_pairs": { + "condition": lambda v: isinstance(v, int) and v > 0, + }, + "num_phases": { + "condition": lambda v: isinstance(v, int) and v >= 2, + }, + "motion_setup": { + "condition": lambda v: True, # only existence check + }, + "rotation_sign": { + "condition": lambda v: v in (-1, 1), + }, + "mechanical_offset_deg": { + "condition": lambda v: isinstance(v, (int, float)) or v == "auto", + }, + "slot_number": { + "condition": lambda v: isinstance(v, (int)), + }, + "layers": { + "condition": lambda v: isinstance(v, (int)), + }, + "coil_pitch": { + "condition": lambda v: isinstance(v, (int)), + }, + "excitation_type": { + "condition": lambda v: v in ("current", "voltage"), + "transform": str.lower, + }, + "phase_base_name": { + "condition": lambda v: True, # only existence check + }, + "phase_index_style": { + "condition": lambda v: v in ("upper", "lower", "number"), + "transform": str.lower, + }, + "custom_phase_names": { + "condition": lambda v: True, # only existence check + }, + "harmonic_orders": { + "transform": lambda v: [v] if isinstance(v, int) else v, + "condition": lambda v: isinstance(v, list) + and len(v) > 0 + and all(isinstance(h, int) and h > 0 for h in v) + }, + "include_inductances": { + "condition": lambda v: True, # only existence check + }, + "inductance_mode": { + "condition": lambda v: v in ("full", "main"), + "transform": str.lower, + }, + "resistance_stat": { + "condition": lambda v: isinstance(v, (int, float)), + }, + "inductance_endw": { + "condition": lambda v: isinstance(v, (int, float)), + }, + } + + # ============================================================ + # Initialization + # ============================================================ + + def __init__(self, user_config: dict): + + if not isinstance(user_config, dict): + raise TypeError("Configuration must be a dictionary.") + + self._raw = user_config + self.params = {} + + print("\n--- Checking configuration ---\n") + + self._validate() + + if self.params["mechanical_offset_deg"] == "auto": + self._validate_winding() + + # ============================================================ + # Parameter validation + # ============================================================ + + def _validate(self): + """ + Iterate over all validation rules and validate parameters. + Missing parameters are automatically filled with default values. + """ + for key, rule in self.VALIDATION_RULES.items(): + + condition = rule["condition"] + transform = rule.get("transform") + + self._validate_general(key, condition, transform) + + # Print final configuration + self._print_final_config() + self._build_phase_names() + + def _validate_general(self, key, condition, transform=None): + """ + General validation function. + Steps: + 1. Get parameter value or use default if the parameter is missing + 2. Apply optional transformation (e.g. string normalization) + 3. Check validity condition + 4. If invalid → replace with default value + """ + + # Check if the parameter exists + if key not in self._raw: + value = self.DEFAULTS[key] + print( + f"[INFO] Parameter '{key}' not provided.\n" + f" Using default value: {value}\n" + ) + else: + value=self._raw[key] + + # Apply transformation if defined + if transform is not None: + try: + value = transform(value) + except Exception: + return self._use_default(key, value) + + # Validate the value + if not condition(value): + return self._use_default(key, value) + + # Store validated value + self.params[key] = value + return value + + def _use_default(self, key, value): + """ + Replace invalid parameter value with default. + """ + default = self.DEFAULTS[key] + print( + f"[INFO] Invalid {key}: {value}\n" + f" using default value: {default}\n" + ) + self.params[key] = default + return default + + # ============================================================ + # Phase naming system + # ============================================================ + + def _build_phase_names(self): + m = self.params["num_phases"] + base = self.params["phase_base_name"] + style = self.params["phase_index_style"] + custom = self.params["custom_phase_names"] + + if custom: + self.params["phase_names"] = custom + self.params["phase_codes"] = custom + return + + names = [] + codes = [] + for k in range(m): + if style == "upper": + idx = chr(ord('A') + k) + + elif style == "lower": + idx = chr(ord('a') + k) + + elif style == "number": + idx = str(k + 1) + + else: + raise ValueError("Invalid phase_index_style") + + names.append(f"{base}{idx}") + codes.append(idx) + + self.params["phase_names"] = names + self.params["phase_codes"] = codes + + # ============================================================ + # Winding validation for automatic mechanical offset + # ============================================================ + + def _validate_winding(self): + Q = self.params["slot_number"] + u2 = self.params["layers"] + ys = self.params["coil_pitch"] + + if Q <= 0 or u2 <= 0 or ys <= 0: + print("[INFO] Slot_number, layers and step_winding must be positive integers.\n" + "Automatic calculation of mechanical offset is disabled due to invalid winding parameters. Using 0 deg.") + self.params["mechanical_offset_deg"] = 0.0 + return + + if u2 == 1 and ys % 2 == 0: + print("[INFO] For single-layer winding the step_winding must be odd.\n" + "Automatic calculation of mechanical offset is disabled due to invalid winding parameters. Using 0 deg.") + self.params["mechanical_offset_deg"] = 0.0 + + # ============================================================ + # Configuration summary + # ============================================================ + + def _print_final_config(self): + """ + Print final validated configuration used by the program. + """ + print("\n========== Simulation configuration ==========") + print("\n[INFO] Final configuration used:\n") + for key in (self.params): + print(f" {key}: {self.params[key]}") + print("==============================================\n") \ No newline at end of file diff --git a/Ansys_transforms/write.py b/Ansys_transforms/write.py new file mode 100644 index 0000000..0be9c89 --- /dev/null +++ b/Ansys_transforms/write.py @@ -0,0 +1,68 @@ +def write_aoutvar(variables, inputs, output_dir="."): + """ + Writes variables to an .aoutvar file. + + Parameters + ---------- + variables : list[str] + Lines prepared by builder (builder.variables) + inputs : dict + User inputs containing phases, source and harmonic_orders + output_dir : str + Directory where the file will be written + """ + + m = inputs["num_phases"] + p = inputs["pole_pairs"] + RotSign = inputs["rotation_sign"] + IndMatrix = inputs["inductance_mode"] + if RotSign == -1: + rot = "CW" + else: + rot = "CCW" + + if inputs["excitation_type"] == "voltage": + source = "V" + else: + source = "I" + + harmonics = sorted(inputs["harmonic_orders"]) + + harmonics_str = "".join(str(h) for h in harmonics) + + filename = f"dq_{m}ph_{p}pp_{source}source_h{harmonics_str}_Lmatrix_{IndMatrix}_{rot}.aoutvar" + filepath = f"{output_dir}/{filename}" + + # exports the dictionary to a Python file + filename="output_vars.py" + output_vars = { + name: var["expression"] + for name, var in variables.items() + } + with open(filename, "w") as f: + f.write("output_vars = {\n") + for k, v in output_vars.items(): + f.write(f" '{k}': \"{v}\",\n") + f.write("}\n") + + # create strings for export from dictionary variables to create .aoutvar file + lines = [] + for name, var in variables.items(): + expr = var["expression"] + + if not var["is_constant"]: + expr = f"'{expr}'" + else: + expr = str(expr) + + unit = var["unit"] + db_flag = var["db_flag"] + + line = f"{name} {expr} Double {unit} {db_flag}" + lines.append(line) + + with open(filepath, "w", encoding="utf-8") as f: + for line in lines: + f.write(line + "\n") + + return filepath \ No newline at end of file