From 1d8e60c5df83158122cc1dfea139ac129e82a3f5 Mon Sep 17 00:00:00 2001 From: brentperteet Date: Fri, 13 Feb 2026 11:32:37 -0600 Subject: [PATCH] adding tordoidal simulator and drawing --- README_WEBAPP.md | 83 ++++ app.py | 197 ++++++++++ core.py | 207 ++++++++++ designer.py | 641 +++++++++++++++++++++++++++++++ draw_toroid.py | 368 ++++++++++++++++++ model.py | 26 +- optimizer.py | 7 + requirements.txt | 2 + sim.py | 6 +- sim_toroid.py | 970 +++++++++++++++++++++++++++++++++++++++++++++++ start_webapp.bat | 9 + 11 files changed, 2511 insertions(+), 5 deletions(-) create mode 100644 README_WEBAPP.md create mode 100644 app.py create mode 100644 core.py create mode 100644 designer.py create mode 100644 draw_toroid.py create mode 100644 requirements.txt create mode 100644 sim_toroid.py create mode 100644 start_webapp.bat diff --git a/README_WEBAPP.md b/README_WEBAPP.md new file mode 100644 index 0000000..0792647 --- /dev/null +++ b/README_WEBAPP.md @@ -0,0 +1,83 @@ +# Transformer Optimizer Web App + +A Flask-based web application for optimizing transformer designs and running manual simulations. + +## Pages + +### 1. Optimizer (http://localhost:5000/) +Automatically finds the best tap configuration and input voltage to maximize efficiency. + +### 2. Manual Simulation (http://localhost:5000/manual) +Directly specify tap settings and input voltage to see performance results. + +## Features + +### Optimizer Page +- **Interactive Sliders** for all parameters: + - Load resistance (5Ω to 10kΩ) + - Target power output (1W to 100W) + - Operating frequency (100Hz to 50kHz) + - Core loss (0W to 5W) + - Maximum flux density (0.05T to 0.5T) + - Maximum secondary voltage (10V to 500V) + - Input voltage range (min/max) + +- **Real-time Optimization** using the transformer optimizer + +### Manual Simulation Page +- **Tap Selection**: Choose primary and secondary taps from dropdowns +- **Input Parameters**: Set input voltage and frequency with sliders +- **Load Conditions**: Configure load resistance and core loss + +### Both Pages Display: +- Efficiency (color-coded: green >90%, yellow >80%, red <80%) +- Tap selections with turn counts +- All voltages and currents +- Separate copper losses for primary and secondary +- Flux density and turns ratio +- Input and output power + +## Installation + +1. Install dependencies: +```bash +pip install -r requirements.txt +``` + +## Running the App + +```bash +python app.py +``` + +Then open your browser to: http://localhost:5000 + +## Configuration + +Edit `app.py` and modify the `get_default_transformer()` function to match your transformer specifications: +- Tap positions +- Core area +- Turn counts +- Resistance per turn for each segment + +## API Endpoints + +- `GET /` - Optimizer page +- `GET /manual` - Manual simulation page +- `POST /api/optimize` - Run optimization with parameters +- `POST /api/simulate` - Run manual simulation +- `GET /api/transformer_info` - Get transformer configuration + +## Usage + +### Optimizer Mode +1. Adjust sliders to set your desired operating conditions +2. Click "⚡ Optimize Transformer" +3. View the optimized configuration and performance metrics +4. The optimizer will find the best tap combination and input voltage to maximize efficiency while meeting your constraints + +### Manual Simulation Mode +1. Select primary and secondary taps from dropdowns +2. Adjust input voltage, frequency, load, and core loss with sliders +3. Click "▶️ Run Simulation" +4. View the performance results for your specified configuration diff --git a/app.py b/app.py new file mode 100644 index 0000000..fb3bb8c --- /dev/null +++ b/app.py @@ -0,0 +1,197 @@ +from flask import Flask, render_template, request, jsonify +from model import TransformerModel +from optimizer import TransformerOptimizer + +app = Flask(__name__) + +# Default transformer configuration +# You can modify this to match your actual transformer +def get_default_transformer(): + + primary_taps = [0, 75, 75] + secondary_taps = [0, 100, 150, 150, 150] + + # Resistance per turn for each segment (example values in ohms/turn) + # These would be calculated based on wire gauge, length per turn, etc. + primary_Rp_per_turn = [ + 0.01, + 0.01, + ] + + secondary_Rs_per_turn = [ + 0.004, + 0.024, + 0.024, + 0.024, + ] + + + + tf = TransformerModel( + Ae_mm2=354.0, + Ve_mm3=43900.0, + use_core_loss_model=True, + Np_total=150, + Ns_total=250, + primary_taps=primary_taps, + secondary_taps=secondary_taps, + primary_Rp_per_turn=primary_Rp_per_turn, + secondary_Rs_per_turn=secondary_Rs_per_turn, + ) + return tf + + +@app.route('/') +def index(): + return render_template('index.html') + + +@app.route('/manual') +def manual(): + return render_template('manual.html') + + +@app.route('/api/simulate', methods=['POST']) +def simulate(): + """Run manual simulation with specified tap and voltage""" + try: + data = request.json + + # Extract parameters + primary_tap = int(data.get('primary_tap', 1)) + secondary_tap = int(data.get('secondary_tap', 1)) + Vp_rms = float(data.get('Vp_rms', 12)) + freq_hz = float(data.get('freq_hz', 2000)) + load_ohms = float(data.get('load_ohms', 100)) + + # Create transformer + tf = get_default_transformer() + + # Run simulation (core loss is calculated automatically) + result = tf.simulate( + primary_tap=primary_tap, + secondary_tap=secondary_tap, + Vp_rms=Vp_rms, + freq_hz=freq_hz, + load_ohms=load_ohms, + core_loss_W=0.0, # Will be calculated by model + ) + + return jsonify({ + 'success': True, + 'result': { + 'primary_tap': result['primary_tap'], + 'secondary_tap': result['secondary_tap'], + 'Np_eff': result['Np_eff'], + 'Ns_eff': result['Ns_eff'], + 'Vp_rms': round(result['Vp_rms'], 2), + 'Vs_rms': round(result['Vs_rms'], 2), + 'Ip_rms': round(result['Ip_rms'], 3), + 'Is_rms': round(result['Is_rms'], 3), + 'turns_ratio': round(result['turns_ratio'], 3), + 'P_out_W': round(result['P_out_W'], 2), + 'P_in_W': round(result['P_in_W'], 2), + 'P_cu_W': round(result['P_cu_W'], 2), + 'P_cu_primary_W': round(result['P_cu_primary_W'], 3), + 'P_cu_secondary_W': round(result['P_cu_secondary_W'], 3), + 'P_core_W': round(result['P_core_W'], 2), + 'efficiency': round(result['efficiency'] * 100, 2), + 'B_peak_T': round(result['B_peak_T'], 4), + } + }) + + except Exception as e: + return jsonify({ + 'success': False, + 'error': str(e) + }), 400 + + +@app.route('/api/optimize', methods=['POST']) +def optimize(): + try: + data = request.json + + # Extract parameters + load_ohms = float(data.get('load_ohms', 100)) + target_power_W = float(data.get('target_power_W', 10)) + freq_hz = float(data.get('freq_hz', 2000)) + Vp_min = float(data.get('Vp_min', 5)) + Vp_max = float(data.get('Vp_max', 36)) + Vp_step = float(data.get('Vp_step', 0.5)) + B_max_T = float(data.get('B_max_T', 0.3)) + Vs_max = float(data.get('Vs_max', 200)) + Is_max = float(data.get('Is_max', 1.5)) + power_tolerance_percent = float(data.get('power_tolerance_percent', 2.0)) + + # Create transformer and optimizer + tf = get_default_transformer() + opt = TransformerOptimizer(tf) + + # Run optimization (core loss is calculated automatically) + result = opt.optimize( + load_ohms=load_ohms, + target_power_W=target_power_W, + freq_hz=freq_hz, + Vp_min=Vp_min, + Vp_max=Vp_max, + Vp_step=Vp_step, + B_max_T=B_max_T, + Vs_max=Vs_max, + Is_max=Is_max, + core_loss_W=0.0, # Will be calculated by model + power_tolerance_percent=power_tolerance_percent, + ) + + if result: + return jsonify({ + 'success': True, + 'result': { + 'primary_tap': result.primary_tap, + 'secondary_tap': result.secondary_tap, + 'Np_eff': result.Np_eff, + 'Ns_eff': result.Ns_eff, + 'Vp_rms': round(result.Vp_rms, 2), + 'Vs_rms': round(result.Vs_rms, 2), + 'Ip_rms': round(result.Ip_rms, 3), + 'Is_rms': round(result.Is_rms, 3), + 'turns_ratio': round(result.turns_ratio, 3), + 'P_out_W': round(result.P_out_W, 2), + 'P_in_W': round(result.P_in_W, 2), + 'P_cu_W': round(result.P_cu_W, 2), + 'P_cu_primary_W': round(result.P_cu_primary_W, 3), + 'P_cu_secondary_W': round(result.P_cu_secondary_W, 3), + 'P_core_W': round(result.P_core_W, 2), + 'efficiency': round(result.efficiency * 100, 2), + 'B_peak_T': round(result.B_peak_T, 4), + 'power_error_percent': round(result.power_error_percent, 2), + } + }) + else: + return jsonify({ + 'success': False, + 'error': 'No valid configuration found within constraints' + }) + + except Exception as e: + return jsonify({ + 'success': False, + 'error': str(e) + }), 400 + + +@app.route('/api/transformer_info', methods=['GET']) +def transformer_info(): + """Return transformer configuration information""" + tf = get_default_transformer() + return jsonify({ + 'primary_taps': tf.primary_taps, + 'secondary_taps': tf.secondary_taps, + 'Ae_mm2': tf.Ae_mm2, + 'Np_total': tf.Np_total, + 'Ns_total': tf.Ns_total, + }) + + +if __name__ == '__main__': + app.run(debug=True, host='0.0.0.0', port=5000) diff --git a/core.py b/core.py new file mode 100644 index 0000000..7fd3d91 --- /dev/null +++ b/core.py @@ -0,0 +1,207 @@ +import numpy as np + +# ------------------------------------------------------------------- +# 1. Put your digitized data here +# Units: +# - B in Tesla +# - frequency in kHz (x-axis of your plot) +# - Pv in kW/m^3 (y-axis of your plot) +# +# IMPORTANT: +# The numbers below are ONLY EXAMPLES / PLACEHOLDERS. +# Replace each "f_khz" and "Pv_kW_m3" array with points +# read from your datasheet plot (solid 25 °C lines). +# ------------------------------------------------------------------- + +core_loss_data_25C = { + 0.025: { # 25 mT + "f_khz": np.array([2.93084374470676, + 8.7361542468832, + 44.706096185109, + 125.963436020903, + 193.225112880907, + ], dtype=float), + "Pv_kW_m3": np.array([0.177193262696068, + 0.679670055798598, + 5.03342212801684, + 17.4677120801992, + 30.0788251804309, + ], dtype=float), + }, + 0.050: { # 50 mT + "f_khz": np.array([2.9975926563905, + 9.24208952674728, + 38.618648046052, + 188.92247157063, + ], dtype=float), + "Pv_kW_m3": np.array([0.917768017464649, + 3.47034515205292, + 19.0328184767803, + 123.927503749051, + ], dtype=float), + }, + 0.100: { # 100 mT + "f_khz": np.array([2.48, + 8.84, + 34.12, + 87.85, + 193.23, + ], dtype=float), + "Pv_kW_m3": np.array([3.14, + 13.89, + 69.94, + 216.47, + 564.36, + ], dtype=float), + }, + 0.200: { # 200 mT + "f_khz": np.array([2.14, + 8.64, + 32.25, + 83.05, + 195.41, + ], dtype=float), + "Pv_kW_m3": np.array([10.44, + 56.44, + 276.05, + 891.89, + 2497.56, + ], dtype=float), + }, + 0.300: { # 300 mT + "f_khz": np.array([2.19, + 5.03, + 12.39, + 24.89, + 49.47, + ], dtype=float), + "Pv_kW_m3": np.array([26.83, + 67.97, + 193.07, + 442.55, + 1000, + ], dtype=float), + }, +} + + +# ------------------------------------------------------------------- +# 2. Helper: 1-D log-log interpolation in frequency +# ------------------------------------------------------------------- + +def _loglog_interp(x, xp, fp): + """ + Log-log interpolation: fp vs xp, evaluated at x. + + All arguments are assumed > 0. + """ + x = np.asarray(x, dtype=float) + xp = np.asarray(xp, dtype=float) + fp = np.asarray(fp, dtype=float) + + logx = np.log10(x) + logxp = np.log10(xp) + logfp = np.log10(fp) + + logy = np.interp(logx, logxp, logfp) + return 10.0 ** logy + + +# ------------------------------------------------------------------- +# 3. Main function: Pv(B, f) interpolation +# ------------------------------------------------------------------- + +def core_loss_Pv(B_T, f_Hz, data=core_loss_data_25C): + """ + Interpolate core loss density P_v for given flux density and frequency. + + Parameters + ---------- + B_T : float + Peak flux density in Tesla (e.g. 0.05 for 50 mT). + f_Hz : float or array_like + Frequency in Hz. + data : dict + Dict mapping B_T -> {"f_khz": array, "Pv_kW_m3": array}. + + Returns + ------- + Pv_kW_m3 : float or ndarray + Core loss density in kW/m^3 at 25 °C. + """ + B_T = float(B_T) + f_khz = np.asarray(f_Hz, dtype=float) / 1e3 # plot uses kHz + + if np.any(f_khz <= 0): + raise ValueError("Frequency must be > 0 Hz") + + # Sorted list of available B values from your digitized curves + B_list = np.array(sorted(data.keys())) + + # Handle out-of-range B values + if B_T < B_list.min(): + raise ValueError( + f"B={B_T:.3f} T is below available range " + f"[{B_list.min():.3f}, {B_list.max():.3f}] T" + ) + + # For B > max data (0.3T), extrapolate up to 0.4T using two highest curves + # If B > 0.4T, clamp to the extrapolated value at 0.4T + if B_T > B_list.max(): + # Use the two highest B values for extrapolation + B1 = B_list[-2] # 0.2T + B2 = B_list[-1] # 0.3T + + d1 = data[B1] + d2 = data[B2] + + # Interpolate along frequency (log-log) on each B curve + Pv1 = _loglog_interp(f_khz, d1["f_khz"], d1["Pv_kW_m3"]) + Pv2 = _loglog_interp(f_khz, d2["f_khz"], d2["Pv_kW_m3"]) + + # Extrapolate in log(Pv) vs B (Pv ~ B^b approximately) + logPv1 = np.log10(Pv1) + logPv2 = np.log10(Pv2) + + # Clamp B_T to 0.4T for extrapolation calculation + B_T_clamped = min(B_T, 0.4) + alpha = (B_T_clamped - B1) / (B2 - B1) + logPv = (1.0 - alpha) * logPv1 + alpha * logPv2 + + return 10.0 ** logPv + + # If B_T is (approximately) one of the tabulated curves, just interpolate in f + idx_match = np.where(np.isclose(B_list, B_T, rtol=1e-4, atol=1e-6))[0] + if idx_match.size: + B_key = B_list[idx_match[0]] + d = data[B_key] + return _loglog_interp(f_khz, d["f_khz"], d["Pv_kW_m3"]) + + # Otherwise, find the two B curves that bracket B_T + idx = np.searchsorted(B_list, B_T) + B1 = B_list[idx - 1] + B2 = B_list[idx] + + d1 = data[B1] + d2 = data[B2] + + # Interpolate along frequency (log-log) on each B curve + Pv1 = _loglog_interp(f_khz, d1["f_khz"], d1["Pv_kW_m3"]) + Pv2 = _loglog_interp(f_khz, d2["f_khz"], d2["Pv_kW_m3"]) + + # Now interpolate between B1 and B2. + # Do it in log(Pv) vs B (Pv ~ B^b approximately). + logPv1 = np.log10(Pv1) + logPv2 = np.log10(Pv2) + + alpha = (B_T - B1) / (B2 - B1) + logPv = (1.0 - alpha) * logPv1 + alpha * logPv2 + + return 10.0 ** logPv + + + +if __name__ == '__main__': + p = core_loss_Pv(0.2, 25_000) + print (p) + diff --git a/designer.py b/designer.py new file mode 100644 index 0000000..9cdb2de --- /dev/null +++ b/designer.py @@ -0,0 +1,641 @@ +""" +Toroid Winding Designer +======================= +Determines the feasibility of winding a toroidal transformer given core +geometry, wire gauge, and turn counts. For each winding segment (tap) it +reports: + + - Whether the turns physically fit on the core + - Number of complete layers required and turns per layer + - Total wire length + - DC resistance + - Wire weight + +Wire gauge table +---------------- +Diameters in the AWG table already include insulation (heavy-build magnet +wire), so no additional insulation allowance is applied. + +Geometry assumptions (toroid cross-section is a rectangle w × h) +----------------------------------------------------------------- +Winding builds simultaneously on ALL four faces of the rectangular +cross-section: + - Inner bore radius shrinks by d_wire per layer + - Outer radius grows by d_wire per layer + - Available winding height shrinks by 2*d_wire per layer + +For layer n (0-indexed) the per-turn path length is: + L_turn(n) = π*(OD/2 + n*d + ID/2 - (n+1)*d) + 2*(h - 2*n*d) + = π*((OD + ID)/2 - d) + 2*h - 4*n*d + +Turns that fit in layer n (minimum of inner-bore and height constraints): + turns_inner(n) = floor( π*(ID - 2*(n+1)*d) / d ) + turns_height(n) = floor( (h - 2*n*d) / d ) + turns_per_layer(n) = min(turns_inner(n), turns_height(n)) + +Winding stops when either constraint hits zero (no more room). + +The fill_factor parameter (default 0.35) limits the fraction of the +available window area that can be used. Window area = π*(ID/2)² for a +toroid. Wire cross-section area per turn = π*(d/2)². Max turns from +fill factor = floor(fill_factor * window_area / wire_area). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field +from typing import Optional + +# --------------------------------------------------------------------------- +# Wire gauge table +# AWG: (diameter_mm including insulation, area_mm2 conductor, ohm_per_km) +# --------------------------------------------------------------------------- + +AWG_TABLE: dict[int, tuple[float, float, float]] = { + 20: (0.813, 0.518, 33.292), + 21: (0.724, 0.410, 41.984), + 22: (0.645, 0.326, 52.939), + 23: (0.574, 0.258, 66.781), + 24: (0.511, 0.205, 84.198), + 25: (0.455, 0.162, 106.17), + 26: (0.404, 0.129, 133.86), + 27: (0.361, 0.102, 168.62), + 28: (0.320, 0.081, 212.87), + 29: (0.287, 0.064, 268.4), + 30: (0.254, 0.051, 338.5), +} + +# Copper density in g/mm³ +_COPPER_DENSITY_G_MM3 = 8.96e-3 + +# --------------------------------------------------------------------------- +# Wire cost table +# AWG: (spool_length_m, cost_usd_per_spool, part_number) +# --------------------------------------------------------------------------- + +WIRE_COST_TABLE: dict[int, tuple[float, float, str]] = { + 20: (972, 230.00, "20SNS10"), + 22: (1545, 238.23, "22SNS10"), + 24: (2447, 222.53, "24SNS10"), + 28: (6178, 224.00, "28SNS10"), +} + + +@dataclass +class WireSpec: + """Properties of a single AWG wire.""" + awg: int + diameter_mm: float # including insulation + area_mm2: float # conductor cross-section + ohm_per_km: float # DC resistance + spool_length_m: Optional[float] = None # None if not in cost table + cost_per_spool: Optional[float] = None # USD + part_number: Optional[str] = None + + @property + def radius_mm(self) -> float: + return self.diameter_mm / 2.0 + + @property + def cost_per_m(self) -> Optional[float]: + """USD per metre, or None if cost data unavailable.""" + if self.spool_length_m and self.cost_per_spool: + return self.cost_per_spool / self.spool_length_m + return None + + @classmethod + def from_awg(cls, awg: int) -> "WireSpec": + if awg not in AWG_TABLE: + raise ValueError(f"AWG {awg} not in table. Available: {sorted(AWG_TABLE)}") + d, a, r = AWG_TABLE[awg] + spool_m = cost_usd = pn = None + if awg in WIRE_COST_TABLE: + spool_m, cost_usd, pn = WIRE_COST_TABLE[awg] + return cls(awg=awg, diameter_mm=d, area_mm2=a, ohm_per_km=r, + spool_length_m=spool_m, cost_per_spool=cost_usd, part_number=pn) + + +@dataclass +class ToroidCore: + """ + Toroidal core geometry and magnetic material parameters. + + Parameters + ---------- + ID_mm : float + Inner diameter of the core in mm. + OD_mm : float + Outer diameter of the core in mm. + height_mm : float + Height (axial length) of the core in mm. + Ae_mm2 : float, optional + Effective magnetic cross-section area (mm²). Used for flux-density + calculation. Defaults to cross_section_area_mm2 (geometric value). + Ve_mm3 : float, optional + Effective core volume (mm³). Used for core loss calculation. + Defaults to Ae_mm2 * mean_path_length_mm. Set to 0 to disable + core loss. + Pv_func : callable, optional + Core loss density function with signature:: + + Pv_func(f_hz: float, B_T: float) -> float # kW / m³ + + If None and Ve_mm3 > 0, the simulator falls back to the built-in + ``core.core_loss_Pv`` interpolation (which also returns kW/m³). + Pass ``lambda f, B: 0.0`` to disable core loss while keeping a + non-zero volume. + """ + ID_mm: float + OD_mm: float + height_mm: float + Ae_mm2: Optional[float] = None + Ve_mm3: Optional[float] = None + Pv_func: Optional[object] = None # Callable[[float, float], float] | None + + def __post_init__(self): + if self.ID_mm <= 0 or self.OD_mm <= self.ID_mm or self.height_mm <= 0: + raise ValueError("Require 0 < ID < OD and height > 0") + if self.Ae_mm2 is not None and self.Ae_mm2 <= 0: + raise ValueError("Ae_mm2 must be > 0") + if self.Ve_mm3 is not None and self.Ve_mm3 < 0: + raise ValueError("Ve_mm3 must be >= 0") + + @property + def window_area_mm2(self) -> float: + """Area of the inner bore (the winding window).""" + return math.pi * (self.ID_mm / 2.0) ** 2 + + @property + def cross_section_area_mm2(self) -> float: + """Rectangular cross-section area of the core material.""" + return ((self.OD_mm - self.ID_mm) / 2.0) * self.height_mm + + @property + def mean_path_length_mm(self) -> float: + """Mean magnetic path length through the core centre.""" + return math.pi * (self.OD_mm + self.ID_mm) / 2.0 + + +@dataclass +class WindingSpec: + """ + Specification for one winding (primary or secondary). + + Parameters + ---------- + awg : int or list[int] + Wire gauge(s). A single int applies to all segments. A list must + have one entry per segment (i.e. len(taps) - 1 entries). + Example: awg=26 or awg=[26, 28] for a 2-segment winding. + taps : list[int] + Incremental turns per segment, same convention as model.py. + First element is always 0; each subsequent element is the additional + turns added by that tap segment. + Example: [0, 128, 23] → tap 1 = 128 turns, tap 2 = 128+23 = 151 turns. + name : str + Optional label (e.g. "primary", "secondary"). + """ + awg: "int | list[int]" + taps: list[int] + name: str = "winding" + + def __post_init__(self): + if len(self.taps) < 2: + raise ValueError("taps must have at least 2 elements (first is always 0)") + if self.taps[0] != 0: + raise ValueError("First element of taps must be 0") + if any(t < 0 for t in self.taps): + raise ValueError("All tap increments must be >= 0") + n_segments = len(self.taps) - 1 + if isinstance(self.awg, int): + self.awg = [self.awg] * n_segments + if len(self.awg) != n_segments: + raise ValueError( + f"awg list length ({len(self.awg)}) must match number of segments ({n_segments})" + ) + + @property + def total_turns(self) -> int: + """Total turns for the highest tap.""" + return sum(self.taps) + + def cumulative_turns(self) -> list[int]: + """Cumulative turns at each tap point (1-indexed tap numbers).""" + result = [] + total = 0 + for t in self.taps[1:]: + total += t + result.append(total) + return result + + def segment_turns(self) -> list[int]: + """Turns in each physical segment (same as taps[1:]).""" + return list(self.taps[1:]) + + +@dataclass +class LayerResult: + """Results for one winding layer.""" + layer_index: int # 0-based + turns_capacity: int # max turns that fit in this layer + turns_used: int # turns actually placed in this layer + L_turn_mm: float # mean turn length for this layer (mm) + wire_length_mm: float # total wire length for turns in this layer (mm) + + +@dataclass +class SegmentResult: + """Results for one tap segment.""" + segment_index: int # 0-based + awg: int + turns: int + layers: list[LayerResult] + wire_length_m: float + resistance_ohm: float + weight_g: float + fits: bool # True if all turns fit within fill factor and geometry + + +@dataclass +class WindingResult: + """Full winding analysis result.""" + spec: WindingSpec + core: ToroidCore + fill_factor: float + max_turns_fill: int # turns allowed by fill factor + max_turns_geometry: int # turns allowed by geometry alone + segments: list[SegmentResult] + total_wire_length_m: float + total_resistance_ohm: float + total_weight_g: float + feasible: bool # True if all segments fit + # Cost fields (None when AWG not in cost table) + spools_required: Optional[int] = None + cost_usd: Optional[float] = None + # Internal winding position at end of this winding (for chaining) + _end_layer: int = 0 + _end_turns_in_layer: int = 0 + _end_total_turns: int = 0 + + def summary(self) -> str: + # Header: show AWG as single value if uniform, else show range + awg_list = self.spec.awg # always a list after __post_init__ + awg_str = str(awg_list[0]) if len(set(awg_list)) == 1 else str(awg_list) + lines = [ + f"=== Winding: {self.spec.name} (AWG {awg_str}) ===", + f" Core: ID={self.core.ID_mm}mm OD={self.core.OD_mm}mm H={self.core.height_mm}mm", + f" Fill factor: {self.fill_factor*100:.0f}% " + f"(max turns by fill={self.max_turns_fill}, by geometry={self.max_turns_geometry})", + f" Total turns: {self.spec.total_turns} | Feasible: {'YES' if self.feasible else 'NO -- DOES NOT FIT'}", + "", + ] + cumulative = 0 + for seg in self.segments: + cumulative += seg.turns + tap_num = seg.segment_index + 1 + status = "OK" if seg.fits else "OVERFLOW" + # Show AWG per tap only when mixed + awg_tag = f" AWG {seg.awg}" if len(set(awg_list)) > 1 else "" + lines.append( + f" Tap {tap_num}: {seg.turns:>4} turns " + f"(cumulative {cumulative:>4}){awg_tag} " + f"{seg.wire_length_m:7.2f} m " + f"{seg.resistance_ohm*1000:8.3f} mOhm " + f"{seg.weight_g:7.3f} g [{status}]" + ) + for lr in seg.layers: + lines.append( + f" layer {lr.layer_index:2d}: capacity={lr.turns_capacity:3d} " + f"used={lr.turns_used:3d} " + f"MTL={lr.L_turn_mm:.1f}mm " + f"wire={lr.wire_length_mm/1000:.3f}m" + ) + cost_str = "" + if self.cost_usd is not None: + cost_str = f" | ${self.cost_usd:.2f}" + lines += [ + "", + f" TOTAL: {self.total_wire_length_m:.2f} m | " + f"{self.total_resistance_ohm*1000:.3f} mOhm | " + f"{self.total_weight_g:.2f} g" + + cost_str, + ] + return "\n".join(lines) + + +def _layer_geometry(core: ToroidCore, wire: WireSpec, layer_index: int + ) -> tuple[int, float]: + """ + Calculate the capacity and mean turn length for a given layer index. + + Returns + ------- + (turns_capacity, L_turn_mm) + turns_capacity : maximum turns that fit in this layer (0 if no room) + L_turn_mm : mean path length of one turn in this layer (mm) + """ + n = layer_index + d = wire.diameter_mm + + # Wire centres sit at ID/2 - (n+0.5)*d from the bore axis + wire_centre_radius = core.ID_mm / 2.0 - (n + 0.5) * d + # Available height for this layer (each layer adds one wire diameter + # of build-up on top and bottom faces, reducing the straight section) + avail_height = core.height_mm - 2 * n * d + + if wire_centre_radius <= 0 or avail_height <= 0: + return 0, 0.0 + + # Turns per layer = how many wire diameters fit around the inner + # circumference at the wire centre radius. The height only constrains + # how many *layers* can be wound, not how many turns per layer. + inner_circumference = 2 * math.pi * wire_centre_radius + turns_capacity = int(inner_circumference / d) # floor + + if turns_capacity <= 0: + return 0, 0.0 + + # Mean turn length: inner + outer semicircle + two straight sections + outer_radius = core.OD_mm / 2.0 + (n + 0.5) * d + L_turn_mm = math.pi * (wire_centre_radius + outer_radius) + 2 * avail_height + + return turns_capacity, L_turn_mm + + +def analyse_winding( + core: ToroidCore, + spec: WindingSpec, + fill_factor: float = 0.35, + start_layer: int = 0, + start_turns_in_layer: int = 0, + start_total_turns: int = 0, +) -> WindingResult: + """ + Analyse feasibility and compute wire parameters for a winding. + + Parameters + ---------- + core : ToroidCore + spec : WindingSpec + fill_factor : float + Maximum fraction of the bore window area that may be filled. + Default 0.35 (35%). + start_layer : int + Layer index to begin winding on. Pass the previous winding's + _end_layer to continue mid-layer across windings. + start_turns_in_layer : int + Turns already consumed in start_layer from a previous winding. + start_total_turns : int + Cumulative turns already placed (used for fill-factor accounting + across windings). + + Returns + ------- + WindingResult + """ + # Resolve per-segment wire specs (spec.awg is always a list after __post_init__) + seg_wires = [WireSpec.from_awg(a) for a in spec.awg] + + # For informational max_turns_fill / max_turns_geometry use the thickest + # wire (largest diameter) as a conservative lower bound. + wire_ref = max(seg_wires, key=lambda w: w.diameter_mm) + wire_area_ref = math.pi * (wire_ref.diameter_mm / 2.0) ** 2 + max_turns_fill = int(fill_factor * core.window_area_mm2 / wire_area_ref) + + max_turns_geometry = 0 + layer_idx = 0 + while True: + cap, _ = _layer_geometry(core, wire_ref, layer_idx) + if cap == 0: + break + max_turns_geometry += cap + layer_idx += 1 + + effective_max = max_turns_fill # fill factor is the binding constraint + + segments: list[SegmentResult] = [] + overall_feasible = True + + # Track position across segments: which layer we're on and how many + # turns have been placed in the current layer already. + current_layer = start_layer + turns_in_current_layer = start_turns_in_layer + total_turns_placed = start_total_turns + current_layer_wire_diameter: Optional[float] = None # set on first use + + for seg_idx, (seg_turns, wire) in enumerate(zip(spec.segment_turns(), seg_wires)): + seg_layers: list[LayerResult] = [] + seg_wire_length_mm = 0.0 + seg_fits = True + turns_remaining = seg_turns + + while turns_remaining > 0: + # If the wire gauge changed from what's already on this layer, + # start a fresh layer so gauges are never mixed on one layer. + if (current_layer_wire_diameter is not None + and wire.diameter_mm != current_layer_wire_diameter + and turns_in_current_layer > 0): + current_layer += 1 + turns_in_current_layer = 0 + + cap, L_turn = _layer_geometry(core, wire, current_layer) + current_layer_wire_diameter = wire.diameter_mm + + if cap == 0: + # No more geometry room at all + seg_fits = False + overall_feasible = False + # Record overflow as a single "layer" with 0 capacity + seg_layers.append(LayerResult( + layer_index=current_layer, + turns_capacity=0, + turns_used=turns_remaining, + L_turn_mm=0.0, + wire_length_mm=0.0, + )) + turns_remaining = 0 + break + + available_in_layer = cap - turns_in_current_layer + if available_in_layer <= 0: + # Layer is full, advance + current_layer += 1 + turns_in_current_layer = 0 + continue + + # Check fill-factor cap + if total_turns_placed >= effective_max: + seg_fits = False + overall_feasible = False + seg_layers.append(LayerResult( + layer_index=current_layer, + turns_capacity=available_in_layer, + turns_used=turns_remaining, + L_turn_mm=L_turn, + wire_length_mm=0.0, + )) + turns_remaining = 0 + break + + turns_to_place = min(turns_remaining, + available_in_layer, + effective_max - total_turns_placed) + wire_len = turns_to_place * L_turn + + seg_layers.append(LayerResult( + layer_index=current_layer, + turns_capacity=cap, + turns_used=turns_to_place, + L_turn_mm=L_turn, + wire_length_mm=wire_len, + )) + + seg_wire_length_mm += wire_len + turns_remaining -= turns_to_place + turns_in_current_layer += turns_to_place + total_turns_placed += turns_to_place + + if turns_in_current_layer >= cap: + current_layer += 1 + turns_in_current_layer = 0 + + seg_wire_length_m = seg_wire_length_mm / 1000.0 + seg_resistance = (wire.ohm_per_km / 1000.0) * seg_wire_length_m + # Weight from conductor area × length × density + seg_weight_g = wire.area_mm2 * seg_wire_length_mm * _COPPER_DENSITY_G_MM3 + + segments.append(SegmentResult( + segment_index=seg_idx, + awg=wire.awg, + turns=seg_turns, + layers=seg_layers, + wire_length_m=seg_wire_length_m, + resistance_ohm=seg_resistance, + weight_g=seg_weight_g, + fits=seg_fits, + )) + + total_wire_m = sum(s.wire_length_m for s in segments) + total_resistance = sum(s.resistance_ohm for s in segments) + total_weight = sum(s.weight_g for s in segments) + + # Cost: pro-rated per segment (each may have different AWG/price) + spools_required = None + cost_usd = None + seg_costs = [] + for seg in segments: + w = WireSpec.from_awg(seg.awg) + if w.cost_per_m is not None: + seg_costs.append(w.cost_per_m * seg.wire_length_m) + if seg_costs: + cost_usd = sum(seg_costs) + + return WindingResult( + spec=spec, + core=core, + fill_factor=fill_factor, + max_turns_fill=max_turns_fill, + max_turns_geometry=max_turns_geometry, + segments=segments, + total_wire_length_m=total_wire_m, + total_resistance_ohm=total_resistance, + total_weight_g=total_weight, + feasible=overall_feasible, + spools_required=spools_required, + cost_usd=cost_usd, + _end_layer=current_layer, + _end_turns_in_layer=turns_in_current_layer, + _end_total_turns=total_turns_placed, + ) + + +def design_transformer( + core: ToroidCore, + windings: list[WindingSpec], + fill_factor: float = 0.35, +) -> list[WindingResult]: + """ + Analyse multiple windings on the same core, packing them continuously + mid-layer. Each winding picks up exactly where the previous one left off. + + The fill_factor budget is shared across all windings: once the total turns + placed (regardless of AWG) reaches fill_factor * window_area / wire_area, + no more turns can be added. Because different AWGs have different wire + areas, the budget is tracked in terms of area consumed. + + Parameters + ---------- + core : ToroidCore + windings : list[WindingSpec] + Windings in winding order. + fill_factor : float + Maximum fraction of bore window area consumed by ALL windings combined. + + Returns + ------- + list[WindingResult] (one per winding, in order) + """ + results: list[WindingResult] = [] + cur_layer = 0 + cur_turns_in_layer = 0 + # Track consumed window area (mm²) across all windings to enforce the + # shared fill-factor budget regardless of AWG. + consumed_area_mm2 = 0.0 + budget_area_mm2 = fill_factor * core.window_area_mm2 + + for spec in windings: + # Use the thickest wire in this spec for area-budget accounting + seg_wires = [WireSpec.from_awg(a) for a in spec.awg] + wire_ref = max(seg_wires, key=lambda w: w.diameter_mm) + wire_area_mm2 = math.pi * (wire_ref.diameter_mm / 2.0) ** 2 + # Express the remaining area budget as an equivalent turn count for + # this winding's reference wire gauge, then back-calculate the + # synthetic start offset so that analyse_winding's fill-factor + # headroom check is correct. + turns_already_equivalent = int(consumed_area_mm2 / wire_area_mm2) + + result = analyse_winding( + core=core, + spec=spec, + fill_factor=fill_factor, + start_layer=cur_layer, + start_turns_in_layer=cur_turns_in_layer, + start_total_turns=turns_already_equivalent, + ) + results.append(result) + + # Advance shared layer position + cur_layer = result._end_layer + cur_turns_in_layer = result._end_turns_in_layer + # Update consumed area using actual turns placed in this winding + turns_placed = result._end_total_turns - turns_already_equivalent + consumed_area_mm2 = min(budget_area_mm2, + consumed_area_mm2 + turns_placed * wire_area_mm2) + + return results + + +# --------------------------------------------------------------------------- +# Example / smoke-test +# --------------------------------------------------------------------------- + +if __name__ == "__main__": + core = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8) + + primary = WindingSpec( + awg=[22, 22], # tap 1 = AWG 22, tap 2 = AWG 24 + taps=[0, 25, 50], + name="primary", + ) + + secondary = WindingSpec( + awg=[22, 22, 22, 26], # uniform gauge + taps=[0, 100, 50, 50, 50], + name="secondary", + ) + + results = design_transformer(core, [primary, secondary], fill_factor=0.35) + for result in results: + print(result.summary()) + print() diff --git a/draw_toroid.py b/draw_toroid.py new file mode 100644 index 0000000..9385223 --- /dev/null +++ b/draw_toroid.py @@ -0,0 +1,368 @@ +""" +draw_toroid.py +============== +Top-down (axial) cross-section drawing of a toroidal transformer. + +The view is looking straight down the axis of the toroid. + +What you see in this view +------------------------- + - Grey annulus : the core material (ID to OD) + - White disc : the bore window (diameter = ID) + - Coloured dots : wire cross-sections packed against the inner bore wall, + one concentric ring per winding layer + +Physical basis: each wire turn passes through the bore. Looking down the +axis you see the wire cross-sections arranged in rings inside the bore. +Layer 0 (first layer) sits closest to the bore wall; successive layers sit +inward. Wire centres for layer n are at: + + r(n) = ID/2 - (n + 0.5) * d_wire + +Turns are spaced by d_wire around the circumference of that ring. Up to +turns_capacity dots are drawn per ring (the physical fill of that layer). + +Usage +----- + from designer import ToroidCore, WindingSpec, design_transformer + from draw_toroid import draw_toroid + + core = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8) + primary = WindingSpec(awg=22, taps=[0, 25, 50], name="primary") + sec = WindingSpec(awg=22, taps=[0, 100, 50], name="secondary") + results = design_transformer(core, [primary, secondary]) + + draw_toroid(core, results, "toroid.png") +""" + +from __future__ import annotations + +import math +from typing import List + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.patches as mpatches +from matplotlib.patches import Circle +import numpy as np + +from designer import ToroidCore, WindingResult, WireSpec + + +# --------------------------------------------------------------------------- +# Colour scheme (one per winding) +# --------------------------------------------------------------------------- +_WINDING_COLORS = [ + "#1565C0", # dark blue — winding 0 + "#B71C1C", # dark red — winding 1 + "#2E7D32", # dark green — winding 2 + "#E65100", # dark orange— winding 3 + "#6A1B9A", # dark purple— winding 4 + "#00838F", # dark cyan — winding 5 +] + +_CORE_FACE_COLOR = "#9E9E9E" # core material +_CORE_EDGE_COLOR = "#212121" +_WIRE_EDGE_COLOR = "#000000" +_BG_COLOR = "white" + + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + +def draw_toroid( + core: ToroidCore, + winding_results: List[WindingResult], + output_path: str = "toroid.png", + dpi: int = 180, + fig_size_mm: float = 180.0, +) -> None: + """ + Render a to-scale top-down cross-section of the transformer as a PNG. + + Parameters + ---------- + core : ToroidCore + winding_results : list[WindingResult] + Output of design_transformer() or analyse_winding(). + output_path : str + Destination file (PNG). + dpi : int + Render resolution. + fig_size_mm : float + Square canvas size in mm. + """ + OD = core.OD_mm + ID = core.ID_mm + + # Canvas extent: just enough margin around the OD + margin = OD * 0.12 + extent = OD / 2.0 + margin # mm from centre to canvas edge + + fig_in = fig_size_mm / 25.4 + fig, ax = plt.subplots(figsize=(fig_in, fig_in), facecolor=_BG_COLOR) + ax.set_facecolor(_BG_COLOR) + ax.set_aspect("equal") + ax.set_xlim(-extent, extent) + ax.set_ylim(-extent, extent) + ax.set_xlabel("mm", fontsize=8) + ax.set_ylabel("mm", fontsize=8) + ax.tick_params(labelsize=7) + ax.set_title( + f"Toroid cross-section " + f"ID={ID:.1f} mm OD={OD:.1f} mm H={core.height_mm:.1f} mm", + fontsize=9, pad=6, + ) + + # ----------------------------------------------------------------------- + # Core annulus + # ----------------------------------------------------------------------- + core_ring = mpatches.Annulus( + xy=(0, 0), + r=OD / 2.0, + width=(OD - ID) / 2.0, + facecolor=_CORE_FACE_COLOR, + edgecolor=_CORE_EDGE_COLOR, + linewidth=1.2, + zorder=2, + label=f"Core ID={ID:.1f} OD={OD:.1f} H={core.height_mm:.1f} mm", + ) + ax.add_patch(core_ring) + + # Bore window + bore = Circle( + (0, 0), ID / 2.0, + facecolor=_BG_COLOR, + edgecolor=_CORE_EDGE_COLOR, + linewidth=0.8, + zorder=3, + ) + ax.add_patch(bore) + + + # ----------------------------------------------------------------------- + # Build layer_index -> wire_centre_radius map (accounts for mixed gauges) + # + # Layers are wound in order from the bore wall inward. Each layer + # consumes its own wire diameter of radial space. We walk through all + # layers across all windings in ascending layer_index order and accumulate + # the actual radial depth so the centre radius is correct regardless of + # whether wire gauge changes between segments or windings. + # ----------------------------------------------------------------------- + + # Collect (layer_index, wire_diameter) for every active layer, all windings + _layer_d: dict[int, float] = {} + for wr in winding_results: + for seg in wr.segments: + wire = WireSpec.from_awg(seg.awg) + for lr in seg.layers: + if lr.turns_used > 0 and lr.turns_capacity > 0: + # First winding to touch a layer defines its wire gauge + if lr.layer_index not in _layer_d: + _layer_d[lr.layer_index] = wire.diameter_mm + + # Walk layers in order to accumulate radial depth from the bore wall + layer_centre_r: dict[int, float] = {} + depth = 0.0 + for n in sorted(_layer_d): + d = _layer_d[n] + layer_centre_r[n] = ID / 2.0 - depth - d / 2.0 + depth += d + + # ----------------------------------------------------------------------- + # Pre-compute per-layer start angles (shared across all windings) so + # the total winding arc on each layer is centred at the top (π/2). + # ----------------------------------------------------------------------- + _layer_total_used: dict[int, int] = {} + for _wr in winding_results: + for _seg in _wr.segments: + for _lr in _seg.layers: + if _lr.turns_used > 0 and _lr.turns_capacity > 0: + n_ = _lr.layer_index + _layer_total_used[n_] = _layer_total_used.get(n_, 0) + _lr.turns_used + + # Step = d/r (touching circles), start at 0, wind continuously. + # All segments share the same running angle per layer — they just + # pick up where the previous segment left off and keep going around. + layer_angle_step: dict[int, float] = {} + layer_next_angle: dict[int, float] = {} + for n_ in _layer_total_used: + r_ = layer_centre_r.get(n_, 0.0) + d_ = _layer_d.get(n_, 1.0) + if r_ > 0: + layer_angle_step[n_] = d_ / r_ + layer_next_angle[n_] = 0.0 + + legend_handles: list[mpatches.Patch] = [] + + for w_idx, wr in enumerate(winding_results): + base_color = _WINDING_COLORS[w_idx % len(_WINDING_COLORS)] + awg_list = wr.spec.awg + awg_str = str(awg_list[0]) if len(set(awg_list)) == 1 else f"{awg_list[0]}-{awg_list[-1]}" + n_segs = len(wr.spec.taps) - 1 + feasible_str = "OK" if wr.feasible else "DOES NOT FIT" + legend_handles.append(mpatches.Patch( + facecolor=base_color, + edgecolor=_WIRE_EDGE_COLOR, + linewidth=0.8, + label=f"{wr.spec.name} AWG {awg_str} {n_segs} seg(s) " + f"{wr.spec.total_turns} turns [{feasible_str}]", + )) + + # Colour shading: outermost layers → light, innermost → dark. + for seg in wr.segments: + wire = WireSpec.from_awg(seg.awg) + d = wire.diameter_mm + + # Shade by segment index: seg 0 → lightest, last seg → base colour. + # Use full 0..1 range across all segments for maximum separation. + frac = seg.segment_index / max(1, n_segs - 1) + seg_color = _lighten(base_color, 0.55 * (1.0 - frac)) + + for lr in seg.layers: + if lr.turns_used == 0 or lr.turns_capacity == 0: + continue + + n = lr.layer_index + r = layer_centre_r.get(n, 0.0) + if r <= 0: + continue + + angle_step = layer_angle_step.get(n, 2.0 * math.pi / max(lr.turns_used, 1)) + start_angle = layer_next_angle.get(n, 0.0) + + # Draw at true wire radius. Circles are evenly spaced over + # 360° so they tile the ring; they may overlap on outer + # layers (where arc spacing > wire diameter) or underlap on + # very dense inner layers, but the count and colour are correct. + draw_r = d / 2.0 + + for i in range(lr.turns_used): + a = start_angle + i * angle_step + x = r * math.cos(a) + y = r * math.sin(a) + ax.add_patch(Circle( + (x, y), draw_r, + facecolor=seg_color, + edgecolor=_WIRE_EDGE_COLOR, + linewidth=0.35, + alpha=0.90, + zorder=10 + n, + )) + + # Advance the angle for the next segment on this layer + layer_next_angle[n] = start_angle + lr.turns_used * angle_step + + # Segment legend entry + awg_tag = f" AWG {seg.awg}" if len(set(awg_list)) > 1 else "" + n_active_layers = sum( + 1 for lr in seg.layers if lr.turns_used > 0 and lr.turns_capacity > 0 + ) + legend_handles.append(mpatches.Patch( + facecolor=seg_color, + edgecolor=_WIRE_EDGE_COLOR, + linewidth=0.6, + label=(f" tap {seg.segment_index + 1}: {seg.turns} turns{awg_tag}" + f" {n_active_layers} layer(s)" + f" {seg.resistance_ohm * 1e3:.1f} mOhm"), + )) + + # ----------------------------------------------------------------------- + # Scale bar (bottom-centre) + # ----------------------------------------------------------------------- + bar_mm = _nice_bar(OD) + bx = -bar_mm / 2 + by = -(extent * 0.87) + tick = extent * 0.012 + ax.plot([bx, bx + bar_mm], [by, by], color="black", lw=1.5, zorder=20) + ax.plot([bx, bx], [by - tick, by + tick], color="black", lw=1.5, zorder=20) + ax.plot([bx + bar_mm, bx + bar_mm], [by - tick, by + tick], color="black", lw=1.5, zorder=20) + ax.text(bx + bar_mm / 2, by - tick * 1.5, f"{bar_mm:.0f} mm", + ha="center", va="top", fontsize=7) + + # ----------------------------------------------------------------------- + # Legend (right of axes) + # ----------------------------------------------------------------------- + ax.legend( + handles=legend_handles, + loc="upper left", + bbox_to_anchor=(1.02, 1.0), + borderaxespad=0, + fontsize=7, + frameon=True, + title="Windings & taps", + title_fontsize=8, + ) + + fig.savefig(output_path, dpi=dpi, bbox_inches="tight", facecolor=_BG_COLOR) + plt.close(fig) + print(f"Saved: {output_path}") + + +# --------------------------------------------------------------------------- +# Geometry helpers +# --------------------------------------------------------------------------- + +def _ring_positions(radius: float, n: int, + d_wire: float = 0.0) -> list[tuple[float, float]]: + """ + Place n wire centres uniformly around a full 360° circle. + + n should equal turns_capacity (floor(2πr/d)), so the spacing between + adjacent centres is approximately d_wire — the circles are tangent with + no visible gap. Always covers the full circumference. + """ + if n <= 0 or radius <= 0: + return [] + angles = np.linspace(0.0, 2 * math.pi, n, endpoint=False) + return [(radius * math.cos(a), radius * math.sin(a)) for a in angles] + + + +def _nice_bar(OD_mm: float) -> float: + """Round scale-bar length ≈ 20% of OD.""" + raw = OD_mm * 0.2 + for v in [1, 2, 5, 10, 20, 25, 50]: + if v >= raw: + return float(v) + return float(round(raw / 10) * 10) + + +def _lighten(hex_color: str, amount: float) -> str: + """ + Mix hex_color towards white by `amount` (0=unchanged, 1=white). + Used to shade tap segments within a winding. + """ + amount = max(0.0, min(1.0, amount)) + h = hex_color.lstrip("#") + r, g, b = (int(h[i:i+2], 16) / 255.0 for i in (0, 2, 4)) + r = r + (1 - r) * amount + g = g + (1 - g) * amount + b = b + (1 - b) * amount + return "#{:02X}{:02X}{:02X}".format(int(r * 255), int(g * 255), int(b * 255)) + + +# --------------------------------------------------------------------------- +# Smoke-test +# --------------------------------------------------------------------------- + +if __name__ == "__main__": + from designer import WindingSpec, design_transformer + + core = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8) + + primary = WindingSpec( + awg=[22, 22], + taps=[0, 25, 50], + name="primary", + ) + secondary = WindingSpec( + awg=[22, 22, 22, 26], + taps=[0, 100, 50, 50, 50], + name="secondary", + ) + + results = design_transformer(core, [primary, secondary]) + draw_toroid(core, results, "toroid.png", dpi=180) diff --git a/model.py b/model.py index 1ff2076..c9a3f18 100644 --- a/model.py +++ b/model.py @@ -1,5 +1,6 @@ from dataclasses import dataclass, field from typing import List, Tuple, Optional, Dict +from core import core_loss_Pv @dataclass @@ -18,6 +19,9 @@ class TransformerModel: # Turn counts (total) Np_total: int Ns_total: int + # Core loss model parameters (optional) + Ve_mm3: float = 0.0 # effective core volume in mm^3 (for core loss calculation) + use_core_loss_model: bool = False # if True, calculate core loss from B and f # Tap positions in turns counted from a common reference (e.g. one end of the bobbin) # Must include 0 and total turns if you want to be able to use full winding. primary_taps: List[int] = field(default_factory=lambda: [0]) @@ -173,8 +177,26 @@ class TransformerModel: # Output power (resistive load) P_out = Vs_rms_ideal**2 / load_ohms + # Core loss calculation + if self.use_core_loss_model and self.Ve_mm3 > 0: + # Calculate core loss from B and f using the model + # Convert volume to m^3 + Ve_m3 = self.Ve_mm3 * 1e-9 # mm^3 to m^3 + + # Get core loss density in kW/m^3 + try: + Pv_kW_m3 = core_loss_Pv(B_peak_T, freq_hz) + # Convert to watts + core_loss_calculated = Pv_kW_m3 * Ve_m3 * 1000 # kW/m^3 * m^3 * 1000 = W + except (ValueError, Exception) as e: + # If core loss model fails (e.g., B or f out of range), fall back to provided value + core_loss_calculated = core_loss_W + else: + # Use provided core loss value + core_loss_calculated = core_loss_W + # Total input power (approx) - P_in = P_out + P_cu_total + core_loss_W + P_in = P_out + P_cu_total + core_loss_calculated efficiency = P_out / P_in if P_in > 0 else 1.0 @@ -191,7 +213,7 @@ class TransformerModel: "P_cu_W": P_cu_total, "P_cu_primary_W": P_cu_p, "P_cu_secondary_W": P_cu_s, - "P_core_W": core_loss_W, + "P_core_W": core_loss_calculated, "P_in_W": P_in, "efficiency": efficiency, "primary_tap": primary_tap, diff --git a/optimizer.py b/optimizer.py index 565605f..de6c0fc 100644 --- a/optimizer.py +++ b/optimizer.py @@ -57,6 +57,7 @@ class TransformerOptimizer: Vp_step: float = 0.5, B_max_T: float = 0.3, Vs_max: float = float('inf'), + Is_max: float = float('inf'), core_loss_W: float = 0.0, power_tolerance_percent: float = 2.0, fallback_max_power: bool = True, @@ -71,6 +72,7 @@ class TransformerOptimizer: - Vp_min, Vp_max, Vp_step: Input voltage search range - B_max_T: Maximum allowed flux density (Tesla) - Vs_max: Maximum allowed secondary voltage (V) + - Is_max: Maximum allowed secondary current (A) - core_loss_W: Core loss at this operating point - power_tolerance_percent: Acceptable power delivery error (%) - fallback_max_power: If True and target cannot be met, find max power with best efficiency @@ -115,6 +117,7 @@ class TransformerOptimizer: Vp_step=Vp_step, B_max_T=B_max_T, Vs_max=Vs_max, + Is_max=Is_max, core_loss_W=core_loss_W, power_tolerance_percent=power_tolerance_percent, fallback_max_power=fallback_max_power, @@ -151,6 +154,7 @@ class TransformerOptimizer: Vp_step: float, B_max_T: float, Vs_max: float, + Is_max: float, core_loss_W: float, power_tolerance_percent: float, fallback_max_power: bool, @@ -195,6 +199,9 @@ class TransformerOptimizer: if sim["Vs_rms"] > Vs_max: continue + if sim["Is_rms"] > Is_max: + continue + power_error = abs(sim["P_out_W"] - target_power_W) / target_power_W * 100 result = OptimizationResult( diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..65919f7 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +flask>=2.3.0 +numpy>=1.24.0 diff --git a/sim.py b/sim.py index 3a58ff7..d4cbec6 100644 --- a/sim.py +++ b/sim.py @@ -86,9 +86,9 @@ opt = TransformerOptimizer(tf) # Find optimal configuration for delivering 10W to a 1000 ohm load at 2 kHz result_opt = opt.optimize( - load_ohms=1500.0, - target_power_W=10.0, - freq_hz=45000.0, + load_ohms=10.0, + target_power_W=25.0, + freq_hz=256.0, Vp_min=1.0, Vp_max=36.0, Vp_step=0.5, diff --git a/sim_toroid.py b/sim_toroid.py new file mode 100644 index 0000000..44e62dd --- /dev/null +++ b/sim_toroid.py @@ -0,0 +1,970 @@ +""" +sim_toroid.py +============= +Simulate a toroidal transformer designed with designer.py. + +Given a pair of WindingResult objects (primary / secondary) produced by +designer.py, this module lets you evaluate any operating point and check +it against a set of constraints. + +Typical usage +------------- + from designer import ToroidCore, WindingSpec, design_transformer + from sim_toroid import ToroidSimulator, SimConstraints + + core = ToroidCore( + ID_mm=21.5, OD_mm=46.5, height_mm=22.8, + Ae_mm2=297.0, # effective cross-section (mm²) + Ve_mm3=18500.0, # effective volume (mm³) + Pv_func=my_pv_func, # Pv_func(f_hz, B_T) -> kW/m³ (or None) + ) + primary = WindingSpec(awg=22, taps=[0, 25, 50], name="primary") + secondary = WindingSpec(awg=22, taps=[0, 100, 50], name="secondary") + pri_result, sec_result = design_transformer(core, [primary, secondary]) + + sim = ToroidSimulator(core=core, primary=pri_result, secondary=sec_result) + + constraints = SimConstraints( + B_max_T=0.3, + Vp_max=50.0, + Vs_max=30.0, + Ip_max=5.0, + Is_max=3.0, + P_out_max_W=50.0, + ) + + result = sim.simulate( + Vp_rms=24.0, + freq_hz=50_000.0, + primary_tap=2, + secondary_tap=1, + Z_load=(10.0, 0.0), # 10 Ω purely resistive + constraints=constraints, + ) + print(result) + +Physics model +------------- +The model uses the standard ideal-transformer approximation with +first-order winding resistance correction: + + * Flux density: B_peak = Vp / (4.44 · Np · Ae · f) + * Ideal secondary voltage: Vs_ideal = Vp · (Ns / Np) + * Secondary circuit: Vs_oc drives (Rs + Z_load) series circuit + Is = Vs_ideal / (Z_load + Rs) + Vs = Is · Z_load + * Primary current reflected from secondary plus magnetising (ignored): + Ip = Is · (Ns / Np) [for turns-ratio correction] + Then added primary drop: + Vp_actual ≈ Vp (Vp is the applied voltage, copper drop on primary + reduces flux slightly; accounted by reducing Vp seen + by the core by Ip·Rp) + * Core loss computed via core.Pv_func(f, B) [kW/m³] · core.Ve_mm3 [mm³] +""" + +from __future__ import annotations + +import cmath +import math +from dataclasses import dataclass, field +from typing import Callable, ClassVar, Optional, Tuple + +from designer import ToroidCore, WindingResult + + +# --------------------------------------------------------------------------- +# Result dataclass +# --------------------------------------------------------------------------- + +@dataclass +class SimResult: + """Operating-point results from ToroidSimulator.simulate().""" + + # Inputs (echoed) + Vp_rms: float + freq_hz: float + primary_tap: int + secondary_tap: int + + # Turns + Np_eff: int + Ns_eff: int + turns_ratio: float # Ns / Np + + # Magnetics + B_peak_T: float # peak flux density in core + + # Voltages (RMS magnitudes) + Vs_rms: float # voltage across load + Vp_rms_applied: float # = Vp_rms (the applied primary voltage) + + # Currents (RMS magnitudes) + Ip_rms: float # primary current + Is_rms: float # secondary current + + # Phase angle of load current relative to secondary voltage (radians) + load_phase_rad: float + + # Winding resistances (ohms) + Rp: float # primary winding resistance for active turns + Rs: float # secondary winding resistance for active turns + + # Power breakdown (watts) + P_out_W: float # real power delivered to load + P_cu_primary_W: float # primary copper loss + P_cu_secondary_W: float # secondary copper loss + P_cu_W: float # total copper loss + P_core_W: float # core loss + P_in_W: float # total input power + efficiency: float # P_out / P_in (0–1) + + # Constraint violations (empty list = all OK) + violations: list[str] = field(default_factory=list) + + @property + def feasible(self) -> bool: + return len(self.violations) == 0 + + def __str__(self) -> str: + lines = [ + f"=== Simulation result (tap P{self.primary_tap}/S{self.secondary_tap}, " + f"{self.freq_hz/1e3:.1f} kHz, Vp={self.Vp_rms:.2f} V) ===", + f" Turns : Np={self.Np_eff} Ns={self.Ns_eff} ratio={self.turns_ratio:.4f}", + f" Flux density : B_peak = {self.B_peak_T*1e3:.2f} mT", + f" Winding R : Rp = {self.Rp*1e3:.3f} mOhm Rs = {self.Rs*1e3:.3f} mOhm", + f" Voltages (RMS) : Vp = {self.Vp_rms_applied:.4f} V Vs = {self.Vs_rms:.4f} V", + f" Currents (RMS) : Ip = {self.Ip_rms:.4f} A Is = {self.Is_rms:.4f} A", + f" Load phase : {math.degrees(self.load_phase_rad):.1f} deg", + f" Power : P_out={self.P_out_W:.4f} W P_cu={self.P_cu_W:.4f} W " + f"P_core={self.P_core_W:.4f} W P_in={self.P_in_W:.4f} W", + f" Efficiency : {self.efficiency*100:.2f}%", + ] + if self.violations: + lines.append(f" *** CONSTRAINT VIOLATIONS: {', '.join(self.violations)} ***") + else: + lines.append(" Constraints : all satisfied") + return "\n".join(lines) + + +# --------------------------------------------------------------------------- +# Constraints dataclass +# --------------------------------------------------------------------------- + +@dataclass +class SimConstraints: + """ + Operating limits. Set any limit to float('inf') / None to disable it. + + Parameters + ---------- + B_max_T : float + Maximum peak flux density (Tesla). Default 0.3 T. + Vp_max : float + Maximum primary RMS voltage. + Vs_max : float + Maximum secondary RMS voltage across load. + Ip_max : float + Maximum primary RMS current. + Is_max : float + Maximum secondary RMS current. + P_out_max_W : float + Maximum output power (W). + """ + B_max_T: float = 0.3 + Vp_max: float = float('inf') + Vs_max: float = float('inf') + Ip_max: float = float('inf') + Is_max: float = float('inf') + P_out_max_W: float = float('inf') + + +# --------------------------------------------------------------------------- +# Simulator +# --------------------------------------------------------------------------- + +class ToroidSimulator: + """ + Simulate an operating point of a toroidal transformer designed with + designer.py. + + Parameters + ---------- + core : ToroidCore + Core geometry and magnetic material parameters. Magnetic fields + (Ae_mm2, Ve_mm3, Pv_func) are read directly from the core object; + if omitted they default to the geometric values. + primary : WindingResult + Output of analyse_winding / design_transformer for the primary. + secondary : WindingResult + Output of analyse_winding / design_transformer for the secondary. + """ + + def __init__( + self, + core: ToroidCore, + primary: WindingResult, + secondary: WindingResult, + ) -> None: + self.core = core + self.primary = primary + self.secondary = secondary + + # ------------------------------------------------------------------ + # Internal helpers + # ------------------------------------------------------------------ + + def _effective_turns(self, winding: WindingResult, tap: int) -> int: + """Cumulative turns up to and including tap (1-indexed).""" + n_taps = len(winding.spec.taps) - 1 # number of valid tap numbers + if tap < 1 or tap > n_taps: + raise ValueError( + f"Tap {tap} out of range 1..{n_taps} for winding '{winding.spec.name}'" + ) + return winding.spec.cumulative_turns()[tap - 1] + + def _winding_resistance(self, winding: WindingResult, tap: int) -> float: + """Total DC resistance for turns up to (and including) tap (ohms).""" + n_taps = len(winding.spec.taps) - 1 + if tap < 1 or tap > n_taps: + raise ValueError( + f"Tap {tap} out of range 1..{n_taps} for winding '{winding.spec.name}'" + ) + # Sum resistance across segments 0..tap-1 + return sum(winding.segments[i].resistance_ohm for i in range(tap)) + + def _effective_Ae_mm2(self) -> float: + """Ae_mm2: use core field if set, else geometric cross-section.""" + if self.core.Ae_mm2 is not None: + return self.core.Ae_mm2 + return self.core.cross_section_area_mm2 + + def _effective_Ve_mm3(self) -> float: + """Ve_mm3: use core field if set, else Ae_mm2 * mean_path_length_mm.""" + if self.core.Ve_mm3 is not None: + return self.core.Ve_mm3 + return self._effective_Ae_mm2() * self.core.mean_path_length_mm + + def _core_loss_W(self, B_peak_T: float, freq_hz: float) -> float: + """Compute core loss in watts. Pv_func returns kW/m³; Ve in mm³.""" + Ve_mm3 = self._effective_Ve_mm3() + if Ve_mm3 <= 0: + return 0.0 + + Ve_m3 = Ve_mm3 * 1e-9 # mm³ → m³ + + if self.core.Pv_func is not None: + # User-supplied function returns kW/m³ + try: + Pv_kW_m3 = float(self.core.Pv_func(freq_hz, B_peak_T)) + except Exception: + return 0.0 + return Pv_kW_m3 * 1e3 * Ve_m3 # kW/m³ → W/m³, then × m³ = W + + # Built-in fallback: core_loss_Pv also returns kW/m³ + try: + from core import core_loss_Pv + Pv_kW_m3 = float(core_loss_Pv(B_peak_T, freq_hz)) + return Pv_kW_m3 * 1e3 * Ve_m3 + except Exception: + return 0.0 + + # ------------------------------------------------------------------ + # Main simulate method + # ------------------------------------------------------------------ + + def simulate( + self, + Vp_rms: float, + freq_hz: float, + primary_tap: int, + secondary_tap: int, + Z_load: Tuple[float, float] = (0.0, 0.0), + constraints: Optional[SimConstraints] = None, + ) -> SimResult: + """ + Compute one operating point. + + Parameters + ---------- + Vp_rms : float + Applied primary RMS voltage (V). + freq_hz : float + Excitation frequency (Hz). + primary_tap : int + Primary tap number (1-indexed). + secondary_tap : int + Secondary tap number (1-indexed). + Z_load : (R, X) tuple + Complex load impedance in ohms. R is the resistive part, + X is the reactive part (positive = inductive, negative = capacitive). + For a purely resistive load use (R, 0.0). + For a complex load the simulator computes real output power as + P_out = |Is|² · R_load. + constraints : SimConstraints, optional + Operating limits. If None, no constraint checking is performed. + + Returns + ------- + SimResult + """ + if Vp_rms <= 0: + raise ValueError("Vp_rms must be > 0") + if freq_hz <= 0: + raise ValueError("freq_hz must be > 0") + + R_load, X_load = float(Z_load[0]), float(Z_load[1]) + Z_load_complex = complex(R_load, X_load) + + if abs(Z_load_complex) == 0: + raise ValueError("Z_load magnitude must be > 0 (short circuit not supported)") + + # --- Turns and resistance ------------------------------------------ + Np = self._effective_turns(self.primary, primary_tap) + Ns = self._effective_turns(self.secondary, secondary_tap) + Rp = self._winding_resistance(self.primary, primary_tap) + Rs = self._winding_resistance(self.secondary, secondary_tap) + + if Np == 0 or Ns == 0: + raise ValueError("Effective turns cannot be zero") + + turns_ratio = Ns / Np # a = Ns/Np + + # --- Flux density --------------------------------------------------- + # The primary voltage applied to the winding minus the resistive drop + # sets the flux. We iterate once: compute ideal Ip, correct Vp seen by + # core, then recompute. + # + # Iteration: + # Pass 0 (ideal): B from full Vp + # Pass 1: B from (Vp - Ip_0 * Rp) where Ip_0 is ideal primary current + # + Ae_m2 = self._effective_Ae_mm2() * 1e-6 + + def _compute_op(Vp_core: float): + """Compute operating point given effective core voltage.""" + B = Vp_core / (4.44 * Np * Ae_m2 * freq_hz) + # Ideal open-circuit secondary voltage + Vs_oc = complex(Vp_core * turns_ratio, 0.0) + # Secondary loop: Vs_oc = Is * (Rs + Z_load) + Z_sec_total = complex(Rs, 0.0) + Z_load_complex + Is_phasor = Vs_oc / Z_sec_total + # Voltage across load + Vs_phasor = Is_phasor * Z_load_complex + # Reflect secondary current to primary (ideal transformer) + Ip_reflected = Is_phasor / turns_ratio # phasor, Amperes + return B, Is_phasor, Vs_phasor, Ip_reflected + + # Pass 0: ideal (ignore primary drop for now) + _, Is0, _, Ip0 = _compute_op(Vp_rms) + Ip0_mag = abs(Ip0) + + # Pass 1: correct for primary winding voltage drop + # Primary drop is Ip * Rp (in phase with current; Rp is real) + # Vp_core ≈ Vp_rms - Ip0 * Rp (phasor subtraction) + # For simplicity treat Ip0 as real-valued magnitude for the correction + # (conservative: subtracts in phase with voltage) + Vp_core = max(Vp_rms - Ip0_mag * Rp, 0.0) + + B_peak_T, Is_phasor, Vs_phasor, Ip_phasor = _compute_op(Vp_core) + + Is_rms = abs(Is_phasor) + Vs_rms_out = abs(Vs_phasor) + Ip_rms = abs(Ip_phasor) + + # Load phase angle (angle of Is relative to Vs across load) + load_phase_rad = cmath.phase(Is_phasor) - cmath.phase(Vs_phasor) + + # --- Power ----------------------------------------------------------- + # Real power out = |Is|² · R_load + P_out = Is_rms ** 2 * R_load + P_cu_p = Ip_rms ** 2 * Rp + P_cu_s = Is_rms ** 2 * Rs + P_cu = P_cu_p + P_cu_s + P_core = self._core_loss_W(B_peak_T, freq_hz) + P_in = P_out + P_cu + P_core + efficiency = P_out / P_in if P_in > 0 else 1.0 + + # --- Constraint checking -------------------------------------------- + violations: list[str] = [] + if constraints is not None: + if B_peak_T > constraints.B_max_T: + violations.append( + f"B_peak={B_peak_T*1e3:.1f}mT > B_max={constraints.B_max_T*1e3:.1f}mT" + ) + if Vp_rms > constraints.Vp_max: + violations.append( + f"Vp={Vp_rms:.2f}V > Vp_max={constraints.Vp_max:.2f}V" + ) + if Vs_rms_out > constraints.Vs_max: + violations.append( + f"Vs={Vs_rms_out:.2f}V > Vs_max={constraints.Vs_max:.2f}V" + ) + if Ip_rms > constraints.Ip_max: + violations.append( + f"Ip={Ip_rms:.4f}A > Ip_max={constraints.Ip_max:.4f}A" + ) + if Is_rms > constraints.Is_max: + violations.append( + f"Is={Is_rms:.4f}A > Is_max={constraints.Is_max:.4f}A" + ) + if P_out > constraints.P_out_max_W: + violations.append( + f"P_out={P_out:.2f}W > P_out_max={constraints.P_out_max_W:.2f}W" + ) + + return SimResult( + Vp_rms=Vp_rms, + freq_hz=freq_hz, + primary_tap=primary_tap, + secondary_tap=secondary_tap, + Np_eff=Np, + Ns_eff=Ns, + turns_ratio=turns_ratio, + B_peak_T=B_peak_T, + Vs_rms=Vs_rms_out, + Vp_rms_applied=Vp_rms, + Ip_rms=Ip_rms, + Is_rms=Is_rms, + load_phase_rad=load_phase_rad, + Rp=Rp, + Rs=Rs, + P_out_W=P_out, + P_cu_primary_W=P_cu_p, + P_cu_secondary_W=P_cu_s, + P_cu_W=P_cu, + P_core_W=P_core, + P_in_W=P_in, + efficiency=efficiency, + violations=violations, + ) + + def optimize( + self, + freq_hz: float, + Z_load: Tuple[float, float], + target_power_W: float, + constraints: SimConstraints, + Vp_min: float = 1.0, + Vp_max: float = 50.0, + Vp_steps: int = 100, + power_tol_pct: float = 2.0, + ) -> Optional["OptimizeResult"]: + """ + Find the tap combination and input voltage that best delivers + ``target_power_W`` to the load while satisfying all constraints. + + The search evaluates every (primary_tap, secondary_tap, Vp) triple. + A candidate is *feasible* if it passes all hard limits in + ``constraints`` (B_max_T, Vp_max, Vs_max, Ip_max, Is_max). + ``constraints.P_out_max_W`` is treated as a hard upper bound too. + + Among feasible candidates: + + 1. If any candidate delivers power within ``power_tol_pct`` % of + ``target_power_W``, the one with the highest efficiency is returned. + 2. Otherwise (fallback), the candidate with power closest to + ``target_power_W`` is returned; ties broken by highest efficiency. + + Parameters + ---------- + freq_hz : float + Operating frequency (Hz). + Z_load : (R, X) tuple + Complex load impedance (ohms). + target_power_W : float + Desired output power (W). + constraints : SimConstraints + Hard operating limits. All fields are enforced. + Vp_min, Vp_max : float + Input voltage search range (V). + Vp_steps : int + Number of voltage steps (linearly spaced, inclusive of endpoints). + power_tol_pct : float + Acceptable power delivery error as % of target. Default 2 %. + + Returns + ------- + OptimizeResult, or None if no feasible point exists at all. + """ + import numpy as np + + if target_power_W <= 0: + raise ValueError("target_power_W must be > 0") + if Vp_min >= Vp_max: + raise ValueError("Vp_min must be < Vp_max") + if Vp_steps < 2: + raise ValueError("Vp_steps must be >= 2") + + voltages = np.linspace(Vp_min, Vp_max, Vp_steps) + n_ptaps = len(self.primary.spec.taps) - 1 + n_staps = len(self.secondary.spec.taps) - 1 + + tol_abs = target_power_W * power_tol_pct / 100.0 + + # Best within tolerance + best_in_tol: Optional[SimResult] = None + best_in_tol_eff = -1.0 + + # Best fallback (closest power, then best efficiency) + best_fallback: Optional[SimResult] = None + best_fallback_err = float('inf') + best_fallback_eff = -1.0 + + for p_tap in range(1, n_ptaps + 1): + for s_tap in range(1, n_staps + 1): + for Vp in voltages: + try: + r = self.simulate( + Vp_rms=float(Vp), + freq_hz=freq_hz, + primary_tap=p_tap, + secondary_tap=s_tap, + Z_load=Z_load, + constraints=constraints, + ) + except ValueError: + continue + + # Skip if any hard constraint violated + if not r.feasible: + continue + + err = abs(r.P_out_W - target_power_W) + + if err <= tol_abs: + if r.efficiency > best_in_tol_eff: + best_in_tol_eff = r.efficiency + best_in_tol = r + else: + if (err < best_fallback_err or + (err == best_fallback_err and r.efficiency > best_fallback_eff)): + best_fallback_err = err + best_fallback_eff = r.efficiency + best_fallback = r + + if best_in_tol is not None: + return OptimizeResult(result=best_in_tol, target_power_W=target_power_W, + power_tol_pct=power_tol_pct, met_target=True) + if best_fallback is not None: + return OptimizeResult(result=best_fallback, target_power_W=target_power_W, + power_tol_pct=power_tol_pct, met_target=False) + return None + + +# --------------------------------------------------------------------------- +# Optimize result dataclass +# --------------------------------------------------------------------------- + +@dataclass +class OptimizeResult: + """ + Result of ToroidSimulator.optimize(). + + Attributes + ---------- + result : SimResult + The best operating point found. + target_power_W : float + The requested target output power. + power_tol_pct : float + The tolerance that was used (%). + met_target : bool + True if result.P_out_W is within power_tol_pct % of target_power_W. + False means the optimizer returned the closest feasible point it + could find (fallback mode). + """ + result: SimResult + target_power_W: float + power_tol_pct: float + met_target: bool + + @property + def power_error_pct(self) -> float: + return abs(self.result.P_out_W - self.target_power_W) / self.target_power_W * 100.0 + + def __str__(self) -> str: + status = "TARGET MET" if self.met_target else "FALLBACK (closest feasible)" + lines = [ + f"=== Optimize result [{status}] ===", + f" Target power : {self.target_power_W:.3f} W " + f"(tol={self.power_tol_pct:.1f}%, actual error={self.power_error_pct:.2f}%)", + str(self.result), + ] + return "\n".join(lines) + + +# --------------------------------------------------------------------------- +# Multi-point sweep: optimise over a grid of frequencies × loads +# --------------------------------------------------------------------------- + +@dataclass +class SweepEntry: + """ + One row of the operating-point sweep dataset. + + Fields that come from the optimizer search inputs are prefixed with no + special marker. Fields derived from the best SimResult are included + directly so the object is self-contained and easy to serialise. + """ + # --- Search inputs ------------------------------------------------------- + freq_hz: float + Z_load_R: float # resistive part of load (ohms) + Z_load_X: float # reactive part of load (ohms) + target_power_W: float + + # --- Optimizer outcome --------------------------------------------------- + met_target: bool # True = within tolerance; False = fallback + power_error_pct: float # |P_out - target| / target * 100 + + # --- Best operating point (None fields mean no feasible solution) -------- + primary_tap: Optional[int] + secondary_tap: Optional[int] + Vp_rms: Optional[float] + Np_eff: Optional[int] + Ns_eff: Optional[int] + turns_ratio: Optional[float] + B_peak_T: Optional[float] + Vs_rms: Optional[float] + Ip_rms: Optional[float] + Is_rms: Optional[float] + load_phase_deg: Optional[float] + Rp_ohm: Optional[float] + Rs_ohm: Optional[float] + P_out_W: Optional[float] + P_cu_W: Optional[float] + P_cu_primary_W: Optional[float] + P_cu_secondary_W: Optional[float] + P_core_W: Optional[float] + P_in_W: Optional[float] + efficiency: Optional[float] + + # CSV column order (class-level constant) + _CSV_FIELDS: ClassVar[list[str]] = [ + "freq_hz", "Z_load_R", "Z_load_X", "target_power_W", + "met_target", "power_error_pct", + "primary_tap", "secondary_tap", "Vp_rms", + "Np_eff", "Ns_eff", "turns_ratio", + "B_peak_T", "Vs_rms", "Ip_rms", "Is_rms", "load_phase_deg", + "Rp_ohm", "Rs_ohm", + "P_out_W", "P_cu_W", "P_cu_primary_W", "P_cu_secondary_W", + "P_core_W", "P_in_W", "efficiency", + ] + + def as_dict(self) -> dict: + return {f: getattr(self, f) for f in self._CSV_FIELDS} + + @staticmethod + def _no_solution( + freq_hz: float, + Z_load: Tuple[float, float], + target_power_W: float, + ) -> "SweepEntry": + """Create a row representing an infeasible / no-solution point.""" + return SweepEntry( + freq_hz=freq_hz, Z_load_R=Z_load[0], Z_load_X=Z_load[1], + target_power_W=target_power_W, + met_target=False, power_error_pct=float('nan'), + primary_tap=None, secondary_tap=None, Vp_rms=None, + Np_eff=None, Ns_eff=None, turns_ratio=None, + B_peak_T=None, Vs_rms=None, Ip_rms=None, Is_rms=None, + load_phase_deg=None, Rp_ohm=None, Rs_ohm=None, + P_out_W=None, P_cu_W=None, P_cu_primary_W=None, + P_cu_secondary_W=None, P_core_W=None, P_in_W=None, + efficiency=None, + ) + + @staticmethod + def _from_opt( + freq_hz: float, + Z_load: Tuple[float, float], + target_power_W: float, + opt: "OptimizeResult", + ) -> "SweepEntry": + r = opt.result + return SweepEntry( + freq_hz=freq_hz, Z_load_R=Z_load[0], Z_load_X=Z_load[1], + target_power_W=target_power_W, + met_target=opt.met_target, + power_error_pct=opt.power_error_pct, + primary_tap=r.primary_tap, + secondary_tap=r.secondary_tap, + Vp_rms=r.Vp_rms, + Np_eff=r.Np_eff, + Ns_eff=r.Ns_eff, + turns_ratio=r.turns_ratio, + B_peak_T=r.B_peak_T, + Vs_rms=r.Vs_rms, + Ip_rms=r.Ip_rms, + Is_rms=r.Is_rms, + load_phase_deg=math.degrees(r.load_phase_rad), + Rp_ohm=r.Rp, + Rs_ohm=r.Rs, + P_out_W=r.P_out_W, + P_cu_W=r.P_cu_W, + P_cu_primary_W=r.P_cu_primary_W, + P_cu_secondary_W=r.P_cu_secondary_W, + P_core_W=r.P_core_W, + P_in_W=r.P_in_W, + efficiency=r.efficiency, + ) + + +def sweep_operating_points( + sim: ToroidSimulator, + frequencies: list[float], + loads: list[Tuple[float, float]], + target_power_W: float, + constraints: SimConstraints, + Vp_min: float = 1.0, + Vp_max: float = 50.0, + Vp_steps: int = 100, + power_tol_pct: float = 2.0, +) -> list[SweepEntry]: + """ + Optimise over a grid of frequencies × complex loads. + + For each (freq, load) combination, calls ``sim.optimize()`` to find the + tap combination and input voltage that best delivers ``target_power_W`` + while satisfying all constraints. + + Parameters + ---------- + sim : ToroidSimulator + frequencies : list[float] + Frequencies to sweep (Hz). + loads : list[(R, X)] + Complex load impedances to sweep (ohms). + target_power_W : float + Desired output power for every operating point. + constraints : SimConstraints + Hard limits applied at every point. + Vp_min, Vp_max : float + Input voltage search range (V). + Vp_steps : int + Number of voltage steps per optimize call. + power_tol_pct : float + Acceptable power delivery error (%). + + Returns + ------- + list[SweepEntry] + One entry per (freq, load) pair, in the order + frequencies × loads (freq varies slowest). + """ + entries: list[SweepEntry] = [] + for freq in frequencies: + for Z in loads: + opt = sim.optimize( + freq_hz=freq, + Z_load=Z, + target_power_W=target_power_W, + constraints=constraints, + Vp_min=Vp_min, + Vp_max=Vp_max, + Vp_steps=Vp_steps, + power_tol_pct=power_tol_pct, + ) + if opt is None: + entries.append(SweepEntry._no_solution(freq, Z, target_power_W)) + else: + entries.append(SweepEntry._from_opt(freq, Z, target_power_W, opt)) + return entries + + +def sweep_to_csv(entries: list[SweepEntry], path: str) -> None: + """ + Write a sweep dataset to a CSV file. + + Parameters + ---------- + entries : list[SweepEntry] + Output of ``sweep_operating_points()``. + path : str + Destination file path (created or overwritten). + """ + import csv + if not entries: + return + fields = SweepEntry._CSV_FIELDS + with open(path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=fields) + writer.writeheader() + for e in entries: + writer.writerow(e.as_dict()) + + +# --------------------------------------------------------------------------- +# Convenience: sweep taps at fixed Vp / freq +# --------------------------------------------------------------------------- + +def sweep_taps( + sim: ToroidSimulator, + Vp_rms: float, + freq_hz: float, + Z_load: Tuple[float, float], + constraints: Optional[SimConstraints] = None, +) -> list[SimResult]: + """ + Evaluate all primary × secondary tap combinations for a fixed operating + point. + + Returns a list of SimResult, sorted by descending efficiency (feasible + results first). + """ + n_ptaps = len(sim.primary.spec.taps) - 1 + n_staps = len(sim.secondary.spec.taps) - 1 + results = [] + for p in range(1, n_ptaps + 1): + for s in range(1, n_staps + 1): + try: + r = sim.simulate( + Vp_rms=Vp_rms, + freq_hz=freq_hz, + primary_tap=p, + secondary_tap=s, + Z_load=Z_load, + constraints=constraints, + ) + results.append(r) + except ValueError: + continue + + # Sort: feasible first, then by descending efficiency + results.sort(key=lambda r: (not r.feasible, -r.efficiency)) + return results + + +# --------------------------------------------------------------------------- +# Example / smoke-test +# --------------------------------------------------------------------------- +def pcv(f_hz, B_T): + pcv_ref = 25 # kW/m^3 + a = 1.3 + b = 2.8 + + ret = pcv_ref * (f_hz / 20000)**a * (B_T / 0.2)**b + return ret + + +if __name__ == "__main__": + + + from designer import WindingSpec, design_transformer + + Ae_mm2 = 142.5 + Le_mm = 106.8 + Ve_mm3 = Ae_mm2 * Le_mm + + core = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8, Ae_mm2=Ae_mm2, Ve_mm3=Ve_mm3, Pv_func=pcv) + + primary = WindingSpec( + awg=[22, 22], + taps=[0, 25, 50], + name="primary", + ) + secondary = WindingSpec( + awg=[22, 22, 22, 26], + taps=[0, 100, 50, 50, 50], + name="secondary", + ) + + pri_result, sec_result = design_transformer(core, [primary, secondary]) + + # Ae_mm2 / Ve_mm3 / Pv_func not set on core → simulator uses geometric defaults + sim = ToroidSimulator(core=core, primary=pri_result, secondary=sec_result) + + constraints = SimConstraints( + B_max_T=1.0, + Vp_max=50.0, + Vs_max=90.0, + Ip_max=5.0, + Is_max=2.0, + P_out_max_W=25.0, + ) + + print("=== Single operating point (10 ohm resistive load) ===") + result = sim.simulate( + Vp_rms=12.0, + freq_hz=256.0, + primary_tap=2, + secondary_tap=1, + Z_load=(100.0, 0.0), + constraints=constraints, + ) + print(result) + print() + + print("=== Complex load (8 ohm + 1 mH at 50 kHz -> X ~= 314 ohm) ===") + X = 2 * math.pi * 50_000.0 * 1e-3 + result2 = sim.simulate( + Vp_rms=24.0, + freq_hz=50_000.0, + primary_tap=2, + secondary_tap=1, + Z_load=(8.0, X), + constraints=constraints, + ) + print(result2) + print() + + print("=== Tap sweep ===") + all_results = sweep_taps(sim, 24.0, 50_000.0, (10.0, 0.0), constraints) + for r in all_results: + feasible = "OK" if r.feasible else "VIOL" + print( + f" P{r.primary_tap}/S{r.secondary_tap} " + f"Np={r.Np_eff:3d} Ns={r.Ns_eff:3d} " + f"Vs={r.Vs_rms:.3f}V Is={r.Is_rms:.4f}A " + f"P_out={r.P_out_W:.3f}W eff={r.efficiency*100:.2f}% [{feasible}]" + ) + + print() + print("=== Optimize: find best taps + Vp to deliver ~10 W at 256 Hz ===") + opt = sim.optimize( + freq_hz=256.0, + Z_load=(10.0, 0.0), + target_power_W=25.0, + constraints=constraints, + Vp_min=1.0, + Vp_max=50.0, + Vp_steps=200, + power_tol_pct=2.0, + ) + if opt is None: + print(" No feasible solution found.") + else: + print(opt) + + print() + print("=== Operating-point sweep (3 freqs x 3 loads) ===") + freqs = [256.0, 870.0, 3140.0, 8900.0] + loads = [(50.0, 0.0), (100.0, 0.0), (200.0, 0.0), (600.0, 0.0)] + entries = sweep_operating_points( + sim=sim, + frequencies=freqs, + loads=loads, + target_power_W=25.0, + constraints=constraints, + Vp_min=1.0, + Vp_max=50.0, + Vp_steps=100, + power_tol_pct=2.0, + ) + # Print a compact table + print(f" {'freq_hz':>8} {'R':>6} {'X':>6} {'met':>5} " + f"{'P_out':>7} {'err%':>5} {'eff%':>6} {'Vp':>6} {'tap':>5}") + for e in entries: + met = "YES" if e.met_target else "NO" + if e.P_out_W is not None: + print(f" {e.freq_hz:>8.1f} {e.Z_load_R:>6.1f} {e.Z_load_X:>6.1f} " + f"{met:>5} {e.P_out_W:>7.3f} {e.power_error_pct:>5.2f} " + f"{e.efficiency*100:>6.2f} {e.Vp_rms:>6.3f} " + f"P{e.primary_tap}/S{e.secondary_tap}") + else: + print(f" {e.freq_hz:>8.1f} {e.Z_load_R:>6.1f} {e.Z_load_X:>6.1f} " + f"{'NO SOL':>5}") + + # Write to CSV + sweep_to_csv(entries, "sweep_results.csv") + print() + print(" CSV written to sweep_results.csv") diff --git a/start_webapp.bat b/start_webapp.bat new file mode 100644 index 0000000..6ee6733 --- /dev/null +++ b/start_webapp.bat @@ -0,0 +1,9 @@ +@echo off +echo Starting Transformer Optimizer Web App... +echo. +echo Open your browser to: http://localhost:5000 +echo. +echo Press Ctrl+C to stop the server +echo. +python app.py +pause