diff --git a/src/gausse/gpu/batch_integrator.py b/src/gausse/gpu/batch_integrator.py index 3c8baac..17e6133 100644 --- a/src/gausse/gpu/batch_integrator.py +++ b/src/gausse/gpu/batch_integrator.py @@ -243,22 +243,25 @@ def integrate_batch_discharge( peak_current = xp.maximum(peak_current, xp.where(active, xp.abs(i_new), 0.0)) past_peak = past_peak | (active & (i_old > 1.0) & (i_new < i_old)) - # ОДИН импульс: обрыв на нуле тока ИЛИ на первом локальном минимуме - # (снаряд начал подкачивать ток обратно). Что раньше. + # ОДИН импульс: обрыв на нуле тока, ИЛИ на первом локальном минимуме + # (снаряд начал подкачивать ток обратно), ИЛИ при падении ниже тока + # удержания ключа (передемпфированный хвост; см. physics/circuit.py). crossed = active & (i_old > 0) & (i_new <= 0) - local_min = active & past_peak & (i_new > i_old) & ~crossed - cut = crossed | local_min + decayed = active & past_peak & (i_new < 0.1) & ~crossed # I_hold=0.1А + local_min = active & past_peak & (i_new > i_old) & ~crossed & ~decayed + cut = crossed | local_min | decayed # zero-crossing: интерполяция к I=0; local_min: берём состояние минимума (old) frac = xp.where(crossed, i_old / (i_old - i_new + 1e-30), 0.0) exit_v = xp.where(crossed, _interp(v_old, v_new, frac), exit_v) exit_x = xp.where(crossed, _interp(x_old, x_new, frac), exit_x) exit_q = xp.where(crossed, _interp(q_old, q_new, frac), exit_q) - exit_v = xp.where(local_min, v_old, exit_v) - exit_x = xp.where(local_min, x_old, exit_x) - exit_q = xp.where(local_min, q_old, exit_q) - # остаточная энергия катушки при обрыве на минимуме -> в потери (freewheel) + cut_at_state = local_min | decayed # обрыв в текущем состоянии (не интерп.) + exit_v = xp.where(cut_at_state, v_old, exit_v) + exit_x = xp.where(cut_at_state, x_old, exit_x) + exit_q = xp.where(cut_at_state, q_old, exit_q) + # остаточная энергия катушки при обрыве -> в потери (freewheel) residual = _magnetic_energy(xp, x_old, i_old, params) - energy_diss = energy_diss + xp.where(local_min, residual, 0.0) + energy_diss = energy_diss + xp.where(cut_at_state, residual, 0.0) committed = committed | cut # валидная коммутация # стиффный конфиг «взорвал» fixed-step (inf/nan) -> стоп, НЕ реализуем diff --git a/src/gausse/gpu/batch_sweep.py b/src/gausse/gpu/batch_sweep.py index 3eaf455..02f380d 100644 --- a/src/gausse/gpu/batch_sweep.py +++ b/src/gausse/gpu/batch_sweep.py @@ -73,6 +73,10 @@ def _prepare_stage(stage, projectile, cur_x_global, cur_v, coil_center_global): if cur_v <= 1e-6: return None, "снаряд остановился/пошёл назад" phys = build_stage_physics(stage, projectile) + # вырожденный контур (см. sim/stage.py) — тот же честный быстрый отказ + tau_worst = phys.inductance_model.l_air_h / (phys.r_total_ohm + phys.r_eddy_coeff_ohm + 1e-12) + if tau_worst < 1e-7: + return None, "вырожденная катушка: L/R < 100 нс (вне области модели)" sensor_global = coil_center_global - stage.sensor_to_coil_distance_m if sensor_global < cur_x_global - 1e-9: return None, "датчик позади снаряда (не сработает)" diff --git a/src/gausse/optim/objective.py b/src/gausse/optim/objective.py index cc745d2..7a03f39 100644 --- a/src/gausse/optim/objective.py +++ b/src/gausse/optim/objective.py @@ -21,7 +21,10 @@ from gausse.physics.inductance import air_core_inductance_wheeler, winding_geome from gausse.sim.coilgun import CoilgunResult, run_coilgun from gausse.storage.schema import RunRecord -MODEL_VERSION = "gausse-physics-v2" # v2: скин/близость (Доуэлл), трение+воздух, паспортные импульсные токи ключей +# v3: коэффициент заполнения железа (fill factor) — закрыт эксплойт «тонкий +# снаряд в толстой катушке»; полная формула отражения вихревого контура; +# гард вырожденных катушек. v2: Доуэлл, трение+воздух, импульсные токи ключей. +MODEL_VERSION = "gausse-physics-v3" @dataclass diff --git a/src/gausse/physics/circuit.py b/src/gausse/physics/circuit.py index f7f6e93..2982ab6 100644 --- a/src/gausse/physics/circuit.py +++ b/src/gausse/physics/circuit.py @@ -100,3 +100,28 @@ def current_local_min_event(t: float, state: np.ndarray, inductance_model, param current_local_min_event.terminal = True current_local_min_event.direction = 1.0 + +# Ток удержания ключа: тиристор сам закрывается, когда ток после пика падает +# ниже I_H (типично 30-100 мА; берём консервативно 0.1 А). Без этого события +# передемпфированный разряд (ток -> 0 асимптотически, нуля не пересекает) +# заставляет решатель молотить весь временной бюджет микросекундными шагами — +# одна конфигурация считалась минутами. Обрыв при I float: + nonlocal peak + i = float(state[I]) + if i > peak: + peak = i + if peak < 2.0 * i_hold_a: # импульса ещё толком не было + return 1.0 + return i - i_hold_a + + event.terminal = True + event.direction = -1.0 + return event diff --git a/src/gausse/physics/inductance.py b/src/gausse/physics/inductance.py index 87d9b73..6f36ea6 100644 --- a/src/gausse/physics/inductance.py +++ b/src/gausse/physics/inductance.py @@ -107,10 +107,15 @@ class CoilInductanceModel: smoothing_width_m: float total_turns: int = 0 b_sat_tesla: float = 1e9 # по умолчанию насыщение отключено (для геом. тестов) + # Доля сечения катушки, занятая железом: A_снаряда/A_среднего витка (≤1). + # λ = n·(B_fe·A_fe + B_возд·(A−A_fe)) ⇒ L = L_air·(1+(μ−1)·fill·overlap). + # БЕЗ этого тонкий снаряд в толстой катушке получал усиление как будто + # железо заполняет ВСЁ сечение — оптимизатор эксплуатировал (КПД «82%»). + iron_fill_factor: float = 1.0 @property def l_iron_coeff(self) -> float: - return self.l_air_h * (self.mu_eff - 1.0) + return self.l_air_h * (self.mu_eff - 1.0) * self.iron_fill_factor @property def saturation_current_a(self) -> float: diff --git a/src/gausse/physics/losses.py b/src/gausse/physics/losses.py index 204a64f..0d2c4e4 100644 --- a/src/gausse/physics/losses.py +++ b/src/gausse/physics/losses.py @@ -46,7 +46,13 @@ def eddy_reflected_resistance_ohm( r_eddy_loop = 2 * math.pi * slug_resistivity_ohm_m / slug_length_m if r_eddy_loop <= 0: return 0.0 - return (char_omega_rad_s * mutual) ** 2 / r_eddy_loop + # собственная индуктивность вихревого контура (виток-соленоид длиной снаряда): + # полная формула отражения R = (ωM)²·R₂/(R₂²+(ωL₂)²) вместо резистивно- + # доминированного приближения (ωM)²/R₂ — иначе на высоких ω (вырожденные + # короткие катушки) отражённое сопротивление нефизично улетало в сотни кОм. + l_eddy_loop = MU_0 * a_slug / slug_length_m + om = char_omega_rad_s + return (om * mutual) ** 2 * r_eddy_loop / (r_eddy_loop**2 + (om * l_eddy_loop) ** 2) def hysteresis_energy_j(slug_radius_m: float, slug_length_m: float, n_cycles: float = 1.0) -> float: diff --git a/src/gausse/sim/stage.py b/src/gausse/sim/stage.py index e9dce72..8c7bcd5 100644 --- a/src/gausse/sim/stage.py +++ b/src/gausse/sim/stage.py @@ -189,6 +189,9 @@ def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> Sta coil_length_m=geometry.coil_length_m, char_omega_rad_s=char_omega, ) + # доля сечения катушки, реально занятая железом (см. inductance.py): + # тонкий снаряд в толстой катушке не может усиливать ВЕСЬ поток + fill = min((projectile.diameter_m / 2) ** 2 / geometry.mean_radius_m**2, 1.0) inductance_model = CoilInductanceModel( l_air_h=l_air_h, coil_length_m=geometry.coil_length_m, @@ -197,6 +200,7 @@ def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> Sta smoothing_width_m=wire_od_m, total_turns=geometry.total_turns, b_sat_tesla=projectile.material.b_sat_tesla, + iron_fill_factor=fill, ) mass_kg = projectile.mass_kg frontal_area_m2 = math.pi * (projectile.diameter_m / 2) ** 2 @@ -237,6 +241,21 @@ def run_stage( inductance_model = phys.inductance_model wire_od_m = stage.wire.insulation_od_mm / 1000 + # Вырожденный контур — вне области применимости модели: при постоянной + # времени L/R (с учётом вихревых при полном перекрытии) короче ~100 нс + # (например, катушка в несколько витков длиной пару мм) ОДУ становится + # неинтегрируемо жёсткой (решатель молотит минутами), а сама инженерная + # модель индуктивности/вихрей там уже не имеет смысла. Честный быстрый отказ. + tau_worst_s = inductance_model.l_air_h / (r_total_ohm + r_eddy_coeff + 1e-12) + if tau_worst_s < 1e-7: + return StageResult( + feasible=False, + reason=( + f"вырожденная катушка: постоянная времени контура {tau_worst_s*1e9:.1f} нс " + f"(< 100 нс) — вне области применимости модели" + ), + ) + x_sensor_m = -stage.sensor_to_coil_distance_m if stage.sensor.kind == "inductive": @@ -298,13 +317,18 @@ def run_stage( (0.0, stage.discharge_time_budget_s), [q0, 0.0, x_fire_m, v_fire_mps], args=(inductance_model, circuit_params), - events=(circuit.zero_current_crossing_event, circuit.current_local_min_event), + events=( + circuit.zero_current_crossing_event, + circuit.current_local_min_event, + # ключ закрывается сам при токе ниже удержания (см. circuit.py) + circuit.make_current_decay_event(), + ), rtol=1e-8, atol=1e-11, ) # какое из событий оборвало разряд первым term_candidates = [] - for ev_idx in range(2): + for ev_idx in range(3): if len(discharge_sol.t_events[ev_idx]) > 0: term_candidates.append((discharge_sol.t_events[ev_idx][0], discharge_sol.y_events[ev_idx][0])) if not term_candidates: diff --git a/tests/test_realism_v2.py b/tests/test_realism_v2.py index 17b7521..26e3f88 100644 --- a/tests/test_realism_v2.py +++ b/tests/test_realism_v2.py @@ -164,3 +164,59 @@ def test_friction_force_smooth_at_zero_velocity(): assert abs(retarding_force_n(p, 1e-6) + retarding_force_n(p, -1e-6)) < 1e-12 # вдали от нуля выходит на полное трение assert retarding_force_n(p, 1.0) > 0.9 * p.retard_const_n + + +def test_iron_fill_factor_punishes_thin_slug_in_fat_coil(): + """Регрессия на эксплойт «КПД 82%»: тонкий снаряд в толстой катушке + не может усиливать весь поток — только долю сечения, которую занимает.""" + from gausse.physics.inductance import CoilInductanceModel + + base = dict(l_air_h=1e-3, coil_length_m=0.05, slug_length_m=0.02, + mu_eff=40.0, smoothing_width_m=0.001, total_turns=300) + full = CoilInductanceModel(**base, iron_fill_factor=1.0) + thin = CoilInductanceModel(**base, iron_fill_factor=0.02) # ⌀4мм в ⌀29мм + # сила при том же токе в 50 раз меньше — халявы больше нет + f_full = float(full.force_newtons(-0.03, 100.0)) + f_thin = float(thin.force_newtons(-0.03, 100.0)) + assert f_thin < 0.03 * f_full + + +def test_degenerate_coil_fails_fast(): + """Регрессия на завис 540с: вырожденная катушка бракуется мгновенно.""" + import time + from dataclasses import replace + + stage = make_stage() + degenerate = replace(stage, turns_per_layer=2, layers=2) # 4 витка, ~2мм + t0 = time.time() + res = run_stage(-0.06, 3.0, degenerate, make_projectile()) + assert time.time() - t0 < 5.0, "вырожденный конфиг должен отваливаться быстро" + # либо честно отбракован гардом, либо честно посчитан быстро — оба исхода ок + if not res.feasible: + assert res.reason + + +def test_eddy_reflection_bounded_at_high_frequency(): + """Полная формула отражения (R₂²+(ωL₂)² в знаменателе): на высоких ω + заметно ниже старого резистивно-доминированного приближения, на низких — + сходится к нему. Совсем вырожденные катушки добивает гард L/R в run_stage.""" + import math as _m + + from gausse.physics.constants import MU_0 + from gausse.physics.losses import EDDY_COUPLING_FACTOR, eddy_reflected_resistance_ohm + + kwargs = dict(slug_radius_m=0.002, slug_length_m=0.03, + slug_resistivity_ohm_m=1.6e-7, mu_eff=259.0, + total_turns=5, coil_length_m=0.002) + + def resistive_approx(om): + a = _m.pi * kwargs["slug_radius_m"] ** 2 + m = EDDY_COUPLING_FACTOR * MU_0 * kwargs["mu_eff"] * (5 / 0.002) * a + r2 = 2 * _m.pi * kwargs["slug_resistivity_ohm_m"] / kwargs["slug_length_m"] + return (om * m) ** 2 / r2 + + om_hi, om_lo = 1.25e5, 1e2 + r_hi = eddy_reflected_resistance_ohm(**kwargs, char_omega_rad_s=om_hi) + r_lo = eddy_reflected_resistance_ohm(**kwargs, char_omega_rad_s=om_lo) + assert r_hi < 0.5 * resistive_approx(om_hi) # высокие ω: ограничена + assert abs(r_lo - resistive_approx(om_lo)) < 0.05 * resistive_approx(om_lo) # низкие: сходится