diff --git a/src/gausse/physics/circuit.py b/src/gausse/physics/circuit.py index e45d901..61701df 100644 --- a/src/gausse/physics/circuit.py +++ b/src/gausse/physics/circuit.py @@ -1,14 +1,14 @@ """ОДУ разряда конденсатора через ключ и катушку, связанное с механикой снаряда. Состояние: [Q, I, x, v] — заряд конденсатора, ток в катушке, положение и -скорость снаряда. Индуктивность зависит от x (проход снаряда), поэтому -уравнение на dI/dt выведено из закона сохранения потокосцепления -Ψ = L(x)*I: +скорость снаряда. Индуктивность железа насыщается по току, поэтому и член +dI/dt, и наведённая ЭДС, и сила выводятся из ОДНОГО потокосцепления +λ(x, I) (см. `inductance.CoilInductanceModel`), что сохраняет энергию: V_c = Q / C dQ/dt = -I - dI/dt = (V_c - I*R_total - I*(dL/dx)*v) / L(x) - dv/dt = F(x, I) / m + dI/dt = (V_c - I*R_total - (∂λ/∂x)*v) / (∂λ/∂I) + dv/dt = F(x, I) / m (F из coenergy того же λ) dx/dt = v """ @@ -16,7 +16,6 @@ from dataclasses import dataclass import numpy as np -from gausse.physics.force import force_on_slug_newtons from gausse.physics.inductance import CoilInductanceModel Q, I, X, V = range(4) @@ -27,10 +26,6 @@ class StageCircuitParams: capacitance_f: float r_total_ohm: float mass_kg: float - mu_eff: float - total_turns: int - coil_length_m: float - b_sat_tesla: float def derivatives( @@ -40,13 +35,13 @@ def derivatives( params: StageCircuitParams, ) -> list[float]: q, current, x, v = state - l_x = inductance_model.l_of_x(x) - dl_dx = inductance_model.dl_dx(x) + dlambda_di = inductance_model.dlambda_di(x, current) + dlambda_dx = inductance_model.dlambda_dx(x, current) 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=dl_dx) + d_i = (v_c - current * params.r_total_ohm - dlambda_dx * v) / dlambda_di + force = inductance_model.force_newtons(x, current) d_v = force / params.mass_kg return [d_q, d_i, v, d_v] diff --git a/src/gausse/physics/constants.py b/src/gausse/physics/constants.py index bce4acc..2fba6a8 100644 --- a/src/gausse/physics/constants.py +++ b/src/gausse/physics/constants.py @@ -4,3 +4,14 @@ MU_0 = 4e-7 * 3.141592653589793 # Гн/м, магнитная проницае RESISTIVITY_COPPER_OHM_M = 1.68e-8 # при ~20°C RESISTIVITY_ALUMINUM_OHM_M = 2.82e-8 # при ~20°C + +# Импульсный (surge) ток намного выше паспортного среднего/непрерывного, +# который лежит в базе компонентов. Для одиночного коротко-импульсного +# выстрела coilgun применяем эти множители к max_current_a детали и +# отбраковываем конфигурацию, если пиковый ток разряда их превышает — +# иначе оптимизатор выдаёт "победителей", которые сожгли бы свой ключ. +# Значения — консервативные типовые (тиристор ITSM ~10× IT(AV), MOSFET +# импульсный ~4× непрерывного, IGBT ~3×, электролит surge ~5× ripple). +# Это оценки, а не паспортные surge-рейтинги конкретных деталей. +SWITCH_SURGE_FACTOR = {"SCR": 10.0, "MOSFET": 4.0, "IGBT": 3.0} +CAPACITOR_SURGE_FACTOR = 5.0 diff --git a/src/gausse/physics/inductance.py b/src/gausse/physics/inductance.py index f74854c..896fb4b 100644 --- a/src/gausse/physics/inductance.py +++ b/src/gausse/physics/inductance.py @@ -12,6 +12,8 @@ from dataclasses import dataclass import numpy as np +from gausse.physics.constants import MU_0 + INCH_M = 0.0254 @@ -76,12 +78,26 @@ def _sigmoid(z: np.ndarray) -> np.ndarray: return 0.5 * (1 + np.tanh(z / 2)) +def _ln_cosh(z: np.ndarray) -> np.ndarray: + """Численно устойчивый ln(cosh z) = |z| + ln((1+e^{-2|z|})/2).""" + az = np.abs(z) + return az + np.log1p(np.exp(-2 * az)) - np.log(2.0) + + @dataclass(frozen=True) class CoilInductanceModel: - """L(x) и аналитическая dL/dx для катушки с проходящим ферромагнитным снарядом. + """Потокосцепление λ(x, I) катушки с ферромагнитным снарядом и насыщением. - x — положение центра снаряда относительно центра катушки (м), катушка - считается центрированной в x=0 в своей локальной системе координат. + x — положение центра снаряда относительно центра катушки (м). + + Вклад железа в потокосцепление насыщается по току: + λ(x, I) = L_air·I + L_iron·overlap(x)·g(I), + g(I) = I_sat·tanh(I/I_sat), L_iron = L_air·(μ_eff−1) + где I_sat — ток, при котором поле в железе достигает B_sat. И сила, и + наведённая ЭДС в контуре выводятся из ОДНОГО и того же λ (метод + coenergy), поэтому энергия сохраняется, а поле не улетает за B_sat. + В пределе малого тока (g(I)≈I) всё сводится к линейной модели + F = 0.5·I²·dL/dx. """ l_air_h: float @@ -89,6 +105,20 @@ class CoilInductanceModel: slug_length_m: float mu_eff: float smoothing_width_m: float + total_turns: int = 0 + b_sat_tesla: float = 1e9 # по умолчанию насыщение отключено (для геом. тестов) + + @property + def l_iron_coeff(self) -> float: + return self.l_air_h * (self.mu_eff - 1.0) + + @property + def saturation_current_a(self) -> float: + """Ток, при котором поле в железе достигает B_sat. inf = без насыщения.""" + if self.total_turns <= 0 or self.mu_eff <= 1.0 or self.b_sat_tesla >= 1e8: + return float("inf") + turns_per_m = self.total_turns / self.coil_length_m + return self.b_sat_tesla / (MU_0 * self.mu_eff * turns_per_m) def _edges(self) -> tuple[float, float]: half_span = (self.coil_length_m + self.slug_length_m) / 2 @@ -106,8 +136,45 @@ class CoilInductanceModel: s2 = _sigmoid((e2 - x) / w) return (s1 * s2 / w) * (s2 - s1) + # --- насыщающиеся токовые функции: g, g', G (интеграл g) --- + def _g(self, i: np.ndarray) -> np.ndarray: + i_sat = self.saturation_current_a + if not np.isfinite(i_sat): + return np.asarray(i, dtype=float) + return i_sat * np.tanh(i / i_sat) + + def _g_prime(self, i: np.ndarray) -> np.ndarray: + i_sat = self.saturation_current_a + if not np.isfinite(i_sat): + return np.ones_like(np.asarray(i, dtype=float)) + return 1.0 - np.tanh(i / i_sat) ** 2 # sech^2, устойчиво + + def _g_integral(self, i: np.ndarray) -> np.ndarray: + i_sat = self.saturation_current_a + if not np.isfinite(i_sat): + return 0.5 * np.asarray(i, dtype=float) ** 2 + return i_sat**2 * _ln_cosh(i / i_sat) + + # --- линейные (ненасыщенные) L(x), dL/dx — для геометрических тестов --- def l_of_x(self, x: np.ndarray) -> np.ndarray: return self.l_air_h * (1 + (self.mu_eff - 1) * self.overlap_fraction(x)) def dl_dx(self, x: np.ndarray) -> np.ndarray: - return self.l_air_h * (self.mu_eff - 1) * self.d_overlap_dx(x) + return self.l_iron_coeff * self.d_overlap_dx(x) + + # --- динамика с насыщением (используется контуром) --- + def dlambda_di(self, x: np.ndarray, i: np.ndarray) -> np.ndarray: + """Дифференциальная индуктивность ∂λ/∂I для члена dI/dt.""" + return self.l_air_h + self.l_iron_coeff * self.overlap_fraction(x) * self._g_prime(i) + + def dlambda_dx(self, x: np.ndarray, i: np.ndarray) -> np.ndarray: + """Пространственный член ∂λ/∂x для наведённой ЭДС.""" + return self.l_iron_coeff * self.d_overlap_dx(x) * self._g(i) + + def force_newtons(self, x: np.ndarray, i: np.ndarray) -> np.ndarray: + """Сила из coenergy: F = ∂/∂x ∫₀ᴵ λ dI' = L_iron·overlap'(x)·G(I).""" + return self.l_iron_coeff * self.d_overlap_dx(x) * self._g_integral(i) + + def effective_field_current_a(self, i: float) -> float: + """'Насыщенный' эффективный ток для оценки поля (поле ~ g(I), кап B_sat).""" + return float(self._g(np.asarray(float(i)))) diff --git a/src/gausse/sim/stage.py b/src/gausse/sim/stage.py index d67b06b..fbdccba 100644 --- a/src/gausse/sim/stage.py +++ b/src/gausse/sim/stage.py @@ -28,6 +28,7 @@ from gausse.physics.inductance import ( effective_permeability, winding_geometry, ) +from gausse.physics.constants import SWITCH_SURGE_FACTOR from gausse.physics.force import saturation_scale, solenoid_field_estimate_tesla from gausse.physics.sensors import ( inductive_signal_peak_v, @@ -144,6 +145,8 @@ def run_stage( slug_length_m=projectile.length_m, mu_eff=mu_eff, smoothing_width_m=wire_od_m, + total_turns=geometry.total_turns, + b_sat_tesla=projectile.material.b_sat_tesla, ) x_sensor_m = -stage.sensor_to_coil_distance_m @@ -193,10 +196,6 @@ def run_stage( 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( @@ -231,12 +230,34 @@ def run_stage( # материала снаряда — численные КПД/скорость в этом случае, вероятно, # завышены относительно реального железа. peak_current_a = float(np.max(np.abs(discharge_sol.y[1]))) + # поле оцениваем по НАСЫЩЕННОМУ эффективному току g(I) — теперь оно не + # улетает за B_sat, а упирается в него (см. inductance.py). + effective_field_current_a = inductance_model.effective_field_current_a(peak_current_a) 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 + mu_eff, geometry.total_turns, geometry.coil_length_m, effective_field_current_a ) + # предупреждение теперь означает реальный заход в насыщение (ток > I_sat) + saturation_warning = peak_current_a > inductance_model.saturation_current_a + + # Жёсткая отбраковка по импульсному току: полупроводниковый ключ реально + # выходит из строя за пределами surge-рейтинга, а электролит перегревается. + # Конфигурация, где пиковый ток разряда это превышает, — не пушка, а + # уничтоженная деталь на первом выстреле, поэтому feasible=False. + # Ключ — жёсткий предел: полупроводник реально выходит из строя за surge. + # Батарея конденсаторов — только предупреждение (см. build_detail): + # электролиты переносят высокий импульсный ток, а в базе — паспортный + # непрерывный/ripple рейтинг, поэтому жёстко браковать по нему нечестно. + switch_surge_a = stage.switch.max_current_a * SWITCH_SURGE_FACTOR.get(stage.switch.kind, 4.0) + if peak_current_a > switch_surge_a: + return StageResult( + feasible=False, + reason=( + f"пиковый ток {peak_current_a:.0f} А превышает импульсный предел ключа " + f"{stage.switch.part_number} ({switch_surge_a:.0f} А = {stage.switch.max_current_a:.0f} А × surge)" + ), + t_sensor_s=t_sensor_s, + t_fire_s=t_sensor_s + fire_delay_s, + ) return StageResult( feasible=True, diff --git a/tests/test_circuit_rlc_analytical.py b/tests/test_circuit_rlc_analytical.py index 612b450..fd6361c 100644 --- a/tests/test_circuit_rlc_analytical.py +++ b/tests/test_circuit_rlc_analytical.py @@ -28,15 +28,12 @@ def _make_frozen_model() -> CoilInductanceModel: def _integrate(r_ohm: float, t_span: tuple[float, float], t_eval: np.ndarray): + # mu_eff=1.0 -> вклад железа нулевой, λ = L_air·I, чистый RLC-контур model = _make_frozen_model() params = StageCircuitParams( capacitance_f=C_F, r_total_ohm=r_ohm, mass_kg=1.0, - mu_eff=1.0, - total_turns=100, - coil_length_m=0.05, - b_sat_tesla=1e9, # эффективно отключает клэмп насыщения для этого теста ) state0 = [Q0, 0.0, 0.0, 0.0] return solve_ivp( diff --git a/tests/test_report.py b/tests/test_report.py index b69e167..52939ae 100644 --- a/tests/test_report.py +++ b/tests/test_report.py @@ -17,14 +17,15 @@ DB = ComponentDatabase.load() def _feasible_config() -> CoilgunConfig: wire = next(w for w in DB.wires if w.part_id == "cu-petv2-1.0mm") capacitor = next(c for c in DB.capacitors if c.part_number == "cap-470uf-400v") - switch = next(s for s in DB.switches if s.part_number == "IRFP250PBF") + # мощный ключ (180А, surge ~720А) — конфигурация с реальным пиковым током + switch = next(s for s in DB.switches if s.part_number == "IRFP4468PBF") sensor = next(s for s in DB.sensors if s.part_number == "TCST2103") steel = next(m for m in DB.projectile_materials if m.name == "Ст3 (конструкционная сталь)") projectile = ProjectileConfig(material=steel, diameter_m=0.008, length_m=0.02) stage = StageConfig( wire=wire, capacitor=capacitor, switch=switch, sensor=sensor, tube_od_m=0.01, - turns_per_layer=20, layers=4, sensor_to_coil_distance_m=0.02, charge_voltage_v=350.0, + turns_per_layer=20, layers=4, sensor_to_coil_distance_m=0.02, charge_voltage_v=200.0, ) return CoilgunConfig( stages=[stage, stage], inter_stage_gaps_m=[0.05], projectile=projectile, diff --git a/tests/test_stage_energy_conservation.py b/tests/test_stage_energy_conservation.py index 5b7137c..0c83741 100644 --- a/tests/test_stage_energy_conservation.py +++ b/tests/test_stage_energy_conservation.py @@ -30,8 +30,9 @@ def _wire(resistivity_ohm_m: float = 1.68e-8) -> WireSpec: 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, + part_number="sw1", kind="MOSFET", max_current_a=500.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", ) @@ -87,29 +88,42 @@ def test_kinetic_energy_increases_for_approaching_slug(): result = _run(_wire(), _switch()) assert result.feasible, result.reason assert result.kinetic_energy_delta_j > 0 - # рост кинетической энергии тут должен означать разгон ВПЕРЁД, а не - # отброс назад (см. test_overpowered_discharge_can_fling_slug_backward) + # рост кинетической энергии тут должен означать разгон ВПЕРЁД assert result.exit_v_mps > 0 -def test_overpowered_discharge_can_fling_slug_backward(): - """Честная фиксация реального режима отказа, а не только "успешных" случаев. - - Слишком большая ёмкость (медленный разряд относительно скорости снаряда) - держит ток высоким уже после того, как снаряд проходит центр катушки — - зона x>0 тянет назад (dL/dx<0), и снаряд может выйти с отрицательной - скоростью. Энергобаланс при этом всё ещё точен (КЭ ~ v^2 не зависит от - знака) — это подтверждает, что находка отражает реальную физику плохо - настроенного расстояния/ёмкости, а не ошибку в уравнениях. +def test_energy_conserved_under_strong_discharge(): + """Энергобаланс должен оставаться точным даже при сильном разряде большой + батареей (нелинейное насыщение железа не должно рвать сохранение энергии — + и сила, и ЭДС выведены из одного λ(x, I), см. inductance.py). """ - overpowered_capacitor = CapacitorSpec( + big_capacitor = CapacitorSpec( part_number="c-big", capacitance_uf=1000.0, voltage_v=400.0, esr_ohm=0.05, max_current_a=300.0, price=300.0, source="test", ) - result = _run(_wire(), _switch(), capacitor=overpowered_capacitor) + result = _run(_wire(), _switch(), capacitor=big_capacitor) assert result.feasible, result.reason - assert result.exit_v_mps < 0 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_field_is_capped_near_saturation(): + """Регрессия на баг '13 Тесла': поле не должно улетать за B_sat железа.""" + strong_capacitor = CapacitorSpec( + part_number="c-strong", capacitance_uf=200.0, voltage_v=600.0, esr_ohm=0.05, + max_current_a=500.0, price=400.0, source="test", + ) + high_v_switch = SwitchSpec( + part_number="sw-hi", kind="MOSFET", max_current_a=1000.0, max_voltage_v=800.0, + on_resistance_ohm=0.01, on_voltage_drop_v=None, turn_on_time_ns=50.0, price=50.0, source="test", + ) + stage = StageConfig( + wire=_wire(), capacitor=strong_capacitor, switch=high_v_switch, sensor=OPTICAL_SENSOR, + tube_od_m=0.01, turns_per_layer=40, layers=5, sensor_to_coil_distance_m=0.02, charge_voltage_v=550.0, + ) + result = run_stage(entry_x_m=-0.05, entry_v_mps=5.0, stage=stage, projectile=PROJECTILE) + assert result.feasible, result.reason + # поле упирается в B_sat снаряда (1.8 Тл в фикстуре), а не улетает в десятки Тл + assert result.peak_field_estimate_tesla <= STEEL.b_sat_tesla * 1.05