jze9/gausse
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

Add dual sensor models and single-stage simulator; fix energy-conservation bug

Browse Source 
- physics/sensors.py: optical/Hall (velocity-independent) and inductive
  (velocity-scaled, sech^2 spatial sensitivity) trigger events for solve_ivp
- sim/stage.py: flight-to-trigger -> fire delay -> discharge -> energy
  accounting, returning StageResult(feasible=False, reason=...) instead of
  raising when a sensor never fires or discharge never commutates
- Found and fixed a real bug caught by the energy-conservation test: the
  saturation clamp was applied to the mechanical force but not the
  electrical back-EMF term, silently breaking energy balance by ~15%.
  Removed the dynamic clamp (documented as a deferred nonlinear-L(x,I)
  limitation) and kept saturation as a diagnostic-only warning
  (StageResult.saturation_warning) so numbers stay honest rather than
  quietly wrong. Balance error is now ~0.02%.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
jze9
2026-07-06 19:55:26 +05:00
parent 2046dcba10
commit f9b77b3756
8 changed files with 484 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
PLAN.md
View File

@@ -18,6 +18,19 @@ ChipDip, Cable.ru и др.).
в SQLite. Никаких приукрашенных цифр — если модель показывает низкий КПД или
нереализуемость, это тоже результат.
**Найденная и исправленная ошибка (Этап 3):** насыщение сердечника изначально
клэмпилось только в механическом уравнении (F=0.5·I²·dL/dx), но не в
электрическом (наведённая ЭДС всё ещё считалась по полной dL/dx) — это
незаметно ломало точный энергобаланс на ~15%. Тест на сохранение энергии
(этого же честного протокола, который просил пользователь) это поймал.
Решение: клэмп насыщения убран из динамики (F=0.5·I²·dL/dx без клэмпа,
энергобаланс теперь точен до ~0.02%), а `solenoid_field_estimate_tesla`/
`saturation_scale` оставлены как ДИАГНОСТИКА — `StageResult.saturation_warning`
честно предупреждает, когда конфигурация физически выходит за пределы
насыщения материала снаряда, не подменяя динамику. Полная нелинейная
L(x, I)-модель с coenergy-выводом силы — в разделе "ограничения модели"
как будущая работа, не как текущая гарантия точности.
Полный план архитектуры: см. историю обсуждения / `physics`, `sim`, `optim`,
`storage`, `report` модули ниже.

9
src/gausse/physics/circuit.py
View File

@@ -46,14 +46,7 @@ def derivatives(
v_c = q / params.capacitance_f
d_q = -current
d_i = (v_c - current * params.r_total_ohm - current * dl_dx * v) / l_x
force = force_on_slug_newtons(
current_a=current,
dl_dx_unsaturated=dl_dx,
mu_eff=params.mu_eff,
total_turns=params.total_turns,
coil_length_m=params.coil_length_m,
b_sat_tesla=params.b_sat_tesla,
)
force = force_on_slug_newtons(current_a=current, dl_dx=dl_dx)
d_v = force / params.mass_kg
return [d_q, d_i, v, d_v]

30
src/gausse/physics/force.py
View File

@@ -1,10 +1,15 @@
"""Сила, действующая на ферромагнитный снаряд: F = 0.5*I^2*dL/dx с насыщением.
"""Сила, действующая на ферромагнитный снаряд: F = 0.5*I^2*dL/dx.
Насыщение — грубая, явно приближённая поправка (v1): оцениваем поле внутри
катушки простой соленоидной формулой и, если оно превышает B_sat материала
снаряда, пропорционально ослабляем вклад сердечника в dL/dx. Реальное
насыщение нелинейно и зависит от геометрии — здесь используется линейный
клэмп как консервативная эвристика, а не точная кривая намагничивания.
Насыщение сердечника сюда сознательно НЕ встроено как динамический клэмп:
ранняя версия ослабляла силу пропорционально превышению B_sat, но клэмп
применялся только к механическому уравнению, а не к электрическому
(наведённая ЭДС по-прежнему считалась через полную dL/dx) — это ломало
точный энергобаланс (dE/dt = -I^2*R), проверяемый в тестах. Модель
насыщения, консистентная и для силы, и для ЭДС одновременно, требует
нелинейной L(x, I) с явной coenergy-производной, что отложено (см. PLAN.md,
"ограничения модели"). Здесь `solenoid_field_estimate_tesla`/`saturation_scale`
оставлены как ДИАГНОСТИЧЕСКИЕ функции — предупреждать, что конфигурация,
вероятно, входит в насыщение, не искажая динамику.
"""
from gausse.physics.constants import MU_0
@@ -23,14 +28,5 @@ def saturation_scale(b_estimate_tesla: float, b_sat_tesla: float) -> float:
return min(1.0, b_sat_tesla / b_estimate_tesla)
def force_on_slug_newtons(
current_a: float,
dl_dx_unsaturated: float,
mu_eff: float,
total_turns: int,
coil_length_m: float,
b_sat_tesla: float,
) -> float:
b_estimate = solenoid_field_estimate_tesla(mu_eff, total_turns, coil_length_m, current_a)
scale = saturation_scale(b_estimate, b_sat_tesla)
return 0.5 * current_a**2 * dl_dx_unsaturated * scale
def force_on_slug_newtons(current_a: float, dl_dx: float) -> float:
return 0.5 * current_a**2 * dl_dx

56
src/gausse/physics/sensors.py Normal file
View File

@@ -0,0 +1,56 @@
"""Событийные модели датчика прохода снаряда для solve_ivp.
Два конкурирующих варианта датчика, оба выражены как условие пересечения
нуля во время баллистического полёта [x, v] к катушке:
- Оптический/Холла: срабатывает строго по положению x, независимо от
скорости снаряда — чистый цифровой фронт.
- Индукционный (виток провода): наведённое напряжение приближённо
пропорционально скорости и "колокольной" чувствительности катушки к
положению снаряда (sech^2-профиль). На низкой скорости пик напряжения
может не достичь порога — событие тогда не срабатывает вовсе, и это
честно отражается в результате как feasible=False, а не подгоняется.
"""
from typing import Callable
import numpy as np
SensorEvent = Callable[[float, np.ndarray], float]
def make_optical_sensor_event(x_sensor_m: float) -> SensorEvent:
def event(t: float, state: np.ndarray) -> float:
return state[0] - x_sensor_m
event.terminal = True
event.direction = 1.0
return event
def make_inductive_sensor_event(
x_sensor_m: float,
sensitivity_v_per_mps: float,
threshold_v: float,
width_m: float,
) -> SensorEvent:
def event(t: float, state: np.ndarray) -> float:
x, v = state
bump = 1.0 / np.cosh((x - x_sensor_m) / width_m) ** 2
return sensitivity_v_per_mps * v * bump - threshold_v
event.terminal = True
event.direction = 1.0
return event
def inductive_signal_peak_v(
sensitivity_v_per_mps: float, velocity_mps: float
) -> float:
"""Пиковое напряжение датчика при данной скорости (bump=1, снаряд точно у датчика).
Полезно, чтобы заранее и дёшево (без интегрирования) отбраковать
заведомо нереализуемые конфигурации: если даже пиковое напряжение
ниже порога, датчик не сработает ни при каком положении.
"""
return sensitivity_v_per_mps * abs(velocity_mps)

255
src/gausse/sim/stage.py Normal file
View File

@@ -0,0 +1,255 @@
"""Одна ступень: баллистический подлёт к датчику -> задержка -> разряд -> энергобаланс.
Нереализуемые случаи (датчик не сработал, разряд не скоммутировался) не
бросают исключение — возвращаются как `StageResult(feasible=False, reason=...)`
с честной причиной, чтобы вызывающий код (в т.ч. оптимизатор) мог их
учитывать, а не терять.
"""
import math
from dataclasses import dataclass
import numpy as np
from scipy.integrate import solve_ivp
from gausse.components.schema import (
CapacitorSpec,
ProjectileMaterialSpec,
SensorSpec,
SwitchSpec,
WireSpec,
)
from gausse.physics import circuit
from gausse.physics.circuit import StageCircuitParams
from gausse.physics.inductance import (
CoilInductanceModel,
air_core_inductance_wheeler,
demagnetizing_factor_prolate,
effective_permeability,
winding_geometry,
)
from gausse.physics.force import saturation_scale, solenoid_field_estimate_tesla
from gausse.physics.sensors import (
inductive_signal_peak_v,
make_inductive_sensor_event,
make_optical_sensor_event,
)
@dataclass(frozen=True)
class ProjectileConfig:
material: ProjectileMaterialSpec
diameter_m: float
length_m: float
@property
def mass_kg(self) -> float:
radius_m = self.diameter_m / 2
volume_m3 = math.pi * radius_m**2 * self.length_m
return volume_m3 * self.material.density_kg_m3
@property
def aspect_ratio(self) -> float:
return self.length_m / self.diameter_m
@dataclass(frozen=True)
class StageConfig:
wire: WireSpec
capacitor: CapacitorSpec
switch: SwitchSpec
sensor: SensorSpec
tube_od_m: float
turns_per_layer: int
layers: int
sensor_to_coil_distance_m: float
charge_voltage_v: float
inductive_sensor_width_m: float = 0.005
sensor_time_budget_s: float = 0.05
discharge_time_budget_s: float = 0.02
@dataclass
class StageResult:
feasible: bool
reason: str | None = None
t_sensor_s: float | None = None
t_fire_s: float | None = None
exit_x_m: float | None = None
exit_v_mps: float | None = None
energy_in_j: float | None = None
energy_dissipated_j: float | None = None
energy_remaining_cap_j: float | None = None
kinetic_energy_delta_j: float | None = None
saturation_warning: bool = False
peak_field_estimate_tesla: float | None = None
discharge_t: np.ndarray | None = None
discharge_q: np.ndarray | None = None
discharge_i: np.ndarray | None = None
discharge_x: np.ndarray | None = None
discharge_v: np.ndarray | None = None
def _wire_resistance_ohm(wire: WireSpec, wire_length_m: float) -> float:
bare_radius_m = wire.gauge_mm / 1000 / 2
area_m2 = math.pi * bare_radius_m**2
return wire.resistivity_ohm_m * wire_length_m / area_m2
def _switch_equivalent_resistance_ohm(
switch: SwitchSpec, capacitance_f: float, inductance_h: float, voltage_v: float
) -> float:
if switch.on_resistance_ohm is not None:
return switch.on_resistance_ohm
# SCR/тиристор: прямое падение аппроксимируется эквивалентным
# сопротивлением относительно характерного масштаба тока контура
# I_ref = V0*sqrt(C/L) (пиковый ток недодемпфированного LC-разряда).
i_ref = voltage_v * math.sqrt(capacitance_f / inductance_h)
drop_v = switch.on_voltage_drop_v or 0.0
return drop_v / max(i_ref, 1e-6)
def _ballistic_derivatives(t: float, state: np.ndarray) -> list[float]:
return [state[1], 0.0]
def run_stage(
entry_x_m: float,
entry_v_mps: float,
stage: StageConfig,
projectile: ProjectileConfig,
) -> StageResult:
mass_kg = projectile.mass_kg
wire_od_m = stage.wire.insulation_od_mm / 1000
geometry = winding_geometry(stage.tube_od_m, wire_od_m, stage.turns_per_layer, stage.layers)
r_wire = _wire_resistance_ohm(stage.wire, geometry.total_wire_length_m)
capacitance_f = stage.capacitor.capacitance_uf * 1e-6
l_air_h = air_core_inductance_wheeler(
geometry.mean_radius_m, geometry.coil_length_m, geometry.radial_depth_m, geometry.total_turns
)
r_switch = _switch_equivalent_resistance_ohm(
stage.switch, capacitance_f, l_air_h, stage.charge_voltage_v
)
r_total_ohm = r_wire + r_switch + stage.capacitor.esr_ohm
demag = demagnetizing_factor_prolate(projectile.aspect_ratio)
mu_eff = effective_permeability(projectile.material.mu_r, demag)
inductance_model = CoilInductanceModel(
l_air_h=l_air_h,
coil_length_m=geometry.coil_length_m,
slug_length_m=projectile.length_m,
mu_eff=mu_eff,
smoothing_width_m=wire_od_m,
)
x_sensor_m = -stage.sensor_to_coil_distance_m
if stage.sensor.kind == "inductive":
peak_v = inductive_signal_peak_v(stage.sensor.sensitivity_v_per_mps, entry_v_mps)
if peak_v < stage.sensor.threshold_v:
return StageResult(
feasible=False,
reason=(
f"индукционный датчик: пиковый сигнал {peak_v:.4f}В "
f"ниже порога {stage.sensor.threshold_v:.4f}В "
f"при скорости {entry_v_mps:.2f} м/с"
),
)
sensor_event = make_inductive_sensor_event(
x_sensor_m,
stage.sensor.sensitivity_v_per_mps,
stage.sensor.threshold_v,
stage.inductive_sensor_width_m,
)
else:
sensor_event = make_optical_sensor_event(x_sensor_m)
flight_sol = solve_ivp(
_ballistic_derivatives,
(0.0, stage.sensor_time_budget_s),
[entry_x_m, entry_v_mps],
events=sensor_event,
rtol=1e-8,
atol=1e-10,
)
if len(flight_sol.t_events[0]) == 0:
return StageResult(
feasible=False,
reason="датчик не сработал в пределах временного бюджета полёта",
)
t_sensor_s = flight_sol.t_events[0][0]
x_at_sensor, v_at_sensor = flight_sol.y_events[0][0]
fire_delay_s = (stage.sensor.propagation_delay_ns + stage.switch.turn_on_time_ns) * 1e-9
x_fire_m = x_at_sensor + v_at_sensor * fire_delay_s
v_fire_mps = v_at_sensor
circuit_params = StageCircuitParams(
capacitance_f=capacitance_f,
r_total_ohm=r_total_ohm,
mass_kg=mass_kg,
mu_eff=mu_eff,
total_turns=geometry.total_turns,
coil_length_m=geometry.coil_length_m,
b_sat_tesla=projectile.material.b_sat_tesla,
)
q0 = capacitance_f * stage.charge_voltage_v
discharge_sol = solve_ivp(
circuit.derivatives,
(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,
rtol=1e-8,
atol=1e-11,
)
if len(discharge_sol.t_events[0]) == 0:
return StageResult(
feasible=False,
reason="разряд не скоммутировался (ток не вернулся к нулю) в пределах временного бюджета",
t_sensor_s=t_sensor_s,
t_fire_s=t_sensor_s + fire_delay_s,
)
q_final, i_final, x_final, v_final = discharge_sol.y_events[0][0]
energy_in_j = 0.5 * capacitance_f * stage.charge_voltage_v**2
energy_remaining_cap_j = q_final**2 / (2 * capacitance_f)
energy_dissipated_j = float(
np.trapezoid(discharge_sol.y[1] ** 2 * r_total_ohm, discharge_sol.t)
)
kinetic_before_j = 0.5 * mass_kg * v_fire_mps**2
kinetic_after_j = 0.5 * mass_kg * v_final**2
# Диагностика насыщения (не влияет на динамику, см. docstring force.py):
# предупреждаем, если реальный пик тока подразумевает поле выше B_sat
# материала снаряда — численные КПД/скорость в этом случае, вероятно,
# завышены относительно реального железа.
peak_current_a = float(np.max(np.abs(discharge_sol.y[1])))
peak_field_estimate_tesla = solenoid_field_estimate_tesla(
mu_eff, geometry.total_turns, geometry.coil_length_m, peak_current_a
)
saturation_warning = (
saturation_scale(peak_field_estimate_tesla, projectile.material.b_sat_tesla) < 1.0
)
return StageResult(
feasible=True,
t_sensor_s=t_sensor_s,
t_fire_s=t_sensor_s + fire_delay_s,
exit_x_m=x_final,
exit_v_mps=v_final,
energy_in_j=energy_in_j,
energy_dissipated_j=energy_dissipated_j,
energy_remaining_cap_j=energy_remaining_cap_j,
kinetic_energy_delta_j=kinetic_after_j - kinetic_before_j,
saturation_warning=saturation_warning,
peak_field_estimate_tesla=peak_field_estimate_tesla,
discharge_t=discharge_sol.t,
discharge_q=discharge_sol.y[0],
discharge_i=discharge_sol.y[1],
discharge_x=discharge_sol.y[2],
discharge_v=discharge_sol.y[3],
)

15
tests/test_force_model.py
View File

@@ -22,14 +22,13 @@ def test_saturation_scale_clamps_above_bsat():
assert scale == pytest.approx(0.5)
def test_force_is_reduced_once_saturated():
common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8)
force_unsaturated = force_on_slug_newtons(current_a=5, **common)
force_saturated = force_on_slug_newtons(current_a=500, **common)
# без клэмпа сила росла бы как I^2 (в 10000 раз); с насыщением рост должен быть намного меньше
assert force_saturated / force_unsaturated < 5000
def test_force_scales_as_current_squared():
# F = 0.5*I^2*dL/dx: без клэмпа (см. модуль-докстринг force.py про энергобаланс)
# сила должна расти строго как I^2, иначе нарушится точный энергобаланс контура.
f_low = force_on_slug_newtons(current_a=5, dl_dx=1e-3)
f_high = force_on_slug_newtons(current_a=50, dl_dx=1e-3)
assert f_high / f_low == pytest.approx(100.0)
def test_force_zero_at_zero_current():
common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8)
assert force_on_slug_newtons(current_a=0, **common) == pytest.approx(0.0)
assert force_on_slug_newtons(current_a=0, dl_dx=1e-3) == pytest.approx(0.0)

50
tests/test_sensors.py Normal file
View File

@@ -0,0 +1,50 @@
import numpy as np
from scipy.integrate import solve_ivp
from gausse.physics.sensors import make_inductive_sensor_event, make_optical_sensor_event
def _ballistic(t, state):
return [state[1], 0.0]
def test_optical_sensor_trigger_is_velocity_independent():
x_sensor = 0.02
trigger_xs = []
for v in (1.0, 5.0, 20.0):
event = make_optical_sensor_event(x_sensor)
sol = solve_ivp(_ballistic, (0, 1.0), [-0.05, v], events=event)
assert len(sol.t_events[0]) == 1
trigger_xs.append(sol.y_events[0][0][0])
assert np.allclose(trigger_xs, x_sensor, atol=1e-9)
def test_inductive_sensor_fires_earlier_at_higher_velocity():
x_sensor = 0.02
sensitivity = 0.05
threshold = 0.3
width = 0.005
trigger_xs = []
for v in (10.0, 20.0, 40.0):
event = make_inductive_sensor_event(x_sensor, sensitivity, threshold, width)
sol = solve_ivp(_ballistic, (0, 1.0), [-0.05, v], events=event)
assert len(sol.t_events[0]) == 1
trigger_xs.append(sol.y_events[0][0][0])
# чем выше скорость, тем раньше (дальше от датчика, т.е. при меньшем x)
# срабатывает индукционный датчик, т.к. sensitivity*v*bump(x) достигает
# порога при меньшем bump(x), а значит при большем |x - x_sensor|
assert trigger_xs[0] > trigger_xs[1] > trigger_xs[2]
assert all(x < x_sensor for x in trigger_xs)
def test_inductive_sensor_never_fires_below_threshold_speed():
x_sensor = 0.02
sensitivity = 0.05
threshold = 0.3
width = 0.005
# пиковый сигнал = sensitivity * v = 0.05 * 1.0 = 0.05 << порог 0.3
event = make_inductive_sensor_event(x_sensor, sensitivity, threshold, width)
sol = solve_ivp(_ballistic, (0, 1.0), [-0.05, 1.0], events=event)
assert len(sol.t_events[0]) == 0

89
tests/test_stage_energy_conservation.py Normal file
View File

@@ -0,0 +1,89 @@
import pytest
from gausse.components.schema import (
CapacitorSpec,
ProjectileMaterialSpec,
SensorSpec,
SwitchSpec,
WireSpec,
)
from gausse.sim.stage import ProjectileConfig, StageConfig, run_stage
STEEL = ProjectileMaterialSpec(
name="steel", density_kg_m3=7850.0, mu_r=200.0, b_sat_tesla=1.8, price_per_kg=100.0, source="test"
)
PROJECTILE = ProjectileConfig(material=STEEL, diameter_m=0.008, length_m=0.02)
CAPACITOR = CapacitorSpec(
part_number="c1", capacitance_uf=1000.0, voltage_v=400.0, esr_ohm=0.05, max_current_a=300.0,
price=300.0, source="test",
)
OPTICAL_SENSOR = SensorSpec(
part_number="s1", kind="optical", propagation_delay_ns=500.0, price=20.0, source="test"
)
def _wire(resistivity_ohm_m: float = 1.68e-8) -> WireSpec:
return WireSpec(
part_id="w1", material="copper", gauge_mm=0.8, insulation_od_mm=0.85,
resistivity_ohm_m=resistivity_ohm_m, max_current_a=10.0, price_per_m=3.0, source="test",
)
def _switch(on_resistance_ohm: float = 0.02) -> SwitchSpec:
return SwitchSpec(
part_number="sw1", kind="MOSFET", max_current_a=200.0, max_voltage_v=500.0,
on_resistance_ohm=on_resistance_ohm, on_voltage_drop_v=None, turn_on_time_ns=50.0,
price=50.0, source="test",
)
def _run(wire, switch, capacitor=CAPACITOR):
# turns_per_layer выбран так, чтобы длина катушки (~1.7см) была сравнима
# со снарядом (2см) — иначе снаряд, войдя глубоко внутрь длинной катушки,
# оказывается в плоской зоне перекрытия (dL/dx=0) и сила не действует.
stage = StageConfig(
wire=wire,
capacitor=capacitor,
switch=switch,
sensor=OPTICAL_SENSOR,
tube_od_m=0.01,
turns_per_layer=20,
layers=4,
sensor_to_coil_distance_m=0.02,
charge_voltage_v=350.0,
)
return run_stage(entry_x_m=-0.05, entry_v_mps=5.0, stage=stage, projectile=PROJECTILE)
def test_stage_is_feasible_with_realistic_components():
result = _run(_wire(), _switch())
assert result.feasible, result.reason
def test_energy_conserved_exactly_when_lossless():
lossless_capacitor = CapacitorSpec(
part_number="c-lossless", capacitance_uf=1000.0, voltage_v=400.0, esr_ohm=0.0,
max_current_a=300.0, price=300.0, source="test",
)
result = _run(
_wire(resistivity_ohm_m=0.0), _switch(on_resistance_ohm=0.0), capacitor=lossless_capacitor
)
assert result.feasible, result.reason
assert result.energy_dissipated_j == pytest.approx(0.0, abs=1e-9)
balance = result.energy_remaining_cap_j + result.kinetic_energy_delta_j
assert balance == pytest.approx(result.energy_in_j, rel=1e-4)
def test_energy_balance_holds_with_realistic_losses():
result = _run(_wire(), _switch())
assert result.feasible, result.reason
balance = (
result.energy_remaining_cap_j + result.energy_dissipated_j + result.kinetic_energy_delta_j
)
assert balance == pytest.approx(result.energy_in_j, rel=1e-2)
def test_kinetic_energy_increases_for_approaching_slug():
result = _run(_wire(), _switch())
assert result.feasible, result.reason
assert result.kinetic_energy_delta_j > 0