Add eddy-current losses in the slug -- the missing loss channel

The user caught the evolution reporting 83.9% efficiency, which is unphysical (real coilguns are single-digit %). Root cause: the model's only loss channel was copper resistance; iron had no losses at all, so the reluctance force came "for free" and the optimizer climbed into that corner. Adds eddy-current loss: the solid steel slug acts as a shorted secondary (1-turn transformer), and its reflected resistance R_eddy = (wM)^2/R_e is added to the circuit while the slug is inside the coil (x overlap(x)). Energy now honestly goes to slug heating instead of kinetic. The formula matches the classical solid-cylinder eddy loss (P ~ sigma*w^2*B^2*a^4), so it's physically grounded, not tuned to a target. Wired through schema/JSON (slug resistivity), losses.py, circuit.py, stage.py, and the GPU batch integrator; energy conservation still holds (0.025%). Effect: best efficiency 83.9% -> ~47%, and the distribution is now realistic (median ~0%, most configs single-digit). 47% is still an optimistic ceiling -- it's the optimizer's single best exploit, and the model still omits tube friction, air drag, skin-effect field penetration (skin depth ~0.44mm < slug radius), and timing imperfection. Hysteresis computed but not in the dynamics (negligible, ~1e-4 J vs eddy). 89 tests. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
2026-07-07 15:03:45 +05:00
parent 4f94811b49
commit 4d16825f7d
7 changed files with 109 additions and 9 deletions
--- a/src/gausse/components/data/projectile_materials.json
+++ b/src/gausse/components/data/projectile_materials.json
@@ -5,7 +5,8 @@
    "mu_r": 550.0,
    "b_sat_tesla": 1.9,
    "price_per_kg": 55.0,
-    "source": "Цена REAL: горячекатаный круг Ст3, ros-met.com/metallurg-moskva.ru, ~49-60₽/кг (одна из позиций ~51.4₽/кг за тонну); плотность/μr/B_sat — ОЦЕНКА (стандартные справочные значения для конструкционной стали, не найдено источника с числами конкретно для Ст3)"
+    "resistivity_ohm_m": 1.6e-7,
+    "source": "Цена REAL: горячекатаный круг Ст3, ros-met.com/metallurg-moskva.ru, ~49-60₽/кг (одна из позиций ~51.4₽/кг за тонну); плотность/μr/B_sat/resistivity — ОЦЕНКА (справочные значения для конструкционной стали; удельное сопротивление ~1.6e-7 Ом·м)"
  },
  {
    "name": "Сталь 10 (низкоуглеродистая)",
@@ -13,7 +14,8 @@
    "mu_r": 1500.0,
    "b_sat_tesla": 2.05,
    "price_per_kg": 55.0,
-    "source": "ОЦЕНКА по всем полям: отдельного розничного объявления на 'Сталь 10' не найдено, цена взята той же полосы, что и Ст3 (~50-70₽/кг); μr выше, чем у Ст3, из-за более высокой чистоты низкоуглеродистой стали — справочная оценка, не измерение"
+    "resistivity_ohm_m": 1.4e-7,
+    "source": "ОЦЕНКА по всем полям: отдельного розничного объявления на 'Сталь 10' не найдено, цена той же полосы, что и Ст3 (~50-70₽/кг); μr выше из-за чистоты; удельное сопротивление ~1.4e-7 Ом·м — справочная оценка"
  },
  {
    "name": "Армко-железо (чистое железо)",
@@ -21,6 +23,7 @@
    "mu_r": 3500.0,
    "b_sat_tesla": 2.15,
    "price_per_kg": 160.0,
-    "source": "ОЦЕНКА по всем полям: специализированный товар, не найдено самостоятельной розничной цены; цена оценена как ~3x от Ст3 на основе общего указания '2-4x дороже обычной стали'; μr/B_sat — стандартные справочные значения для хорошо отожжённого чистого железа"
+    "resistivity_ohm_m": 1.0e-7,
+    "source": "ОЦЕНКА по всем полям: специализированный товар, розничной цены не найдено; цена ~3x от Ст3; μr/B_sat — справочные для отожжённого чистого железа; удельное сопротивление ~1.0e-7 Ом·м (чище -> ниже)"
  }
 ]
--- a/src/gausse/components/schema.py
+++ b/src/gausse/components/schema.py
@@ -99,3 +99,6 @@ class ProjectileMaterialSpec:
    b_sat_tesla: float
    price_per_kg: float
    source: str
+    # электрическое удельное сопротивление — для расчёта вихревых потерь в снаряде
+    # (сплошной проводник в импульсном поле). Дефолт ~ конструкционная сталь.
+    resistivity_ohm_m: float = 1.6e-7
--- a/src/gausse/gpu/batch_integrator.py
+++ b/src/gausse/gpu/batch_integrator.py
@@ -34,6 +34,7 @@ class BatchDischargeParams:
    slug_length_m: "any"
    smoothing_width_m: "any"
    i_sat_a: "any"  # ток насыщения; NO_SATURATION_I_SAT где железа нет
+    r_eddy_coeff_ohm: "any"  # вихревое сопротивление снаряда (× overlap)


 def params_from_models(xp, models, circuit_params) -> BatchDischargeParams:
@@ -60,6 +61,7 @@ def params_from_models(xp, models, circuit_params) -> BatchDischargeParams:
        slug_length_m=col([m.slug_length_m for m in models]),
        smoothing_width_m=col([m.smoothing_width_m for m in models]),
        i_sat_a=col(i_sat),
+        r_eddy_coeff_ohm=col([getattr(c, "r_eddy_coeff_ohm", 0.0) for c in circuit_params]),
    )


@@ -86,9 +88,10 @@ def _derivatives(xp, q, i, x, v, p: BatchDischargeParams):
    dlambda_dx = p.l_iron_coeff * d_overlap * g
    force = p.l_iron_coeff * d_overlap * g_integral

+    r_eff = p.r_total_ohm + p.r_eddy_coeff_ohm * overlap  # вихревые потери снаряда × overlap
    v_c = q / p.capacitance_f
    d_q = -i
-    d_i = (v_c - i * p.r_total_ohm - dlambda_dx * v) / dlambda_di
+    d_i = (v_c - i * r_eff - dlambda_dx * v) / dlambda_di
    d_v = force / p.mass_kg
    d_x = v
    return d_q, d_i, d_x, d_v
--- a/src/gausse/physics/circuit.py
+++ b/src/gausse/physics/circuit.py
@@ -26,6 +26,14 @@ class StageCircuitParams:
    capacitance_f: float
    r_total_ohm: float
    mass_kg: float
+    # отражённое сопротивление вихревых токов снаряда; действует × overlap(x)
+    # (только пока снаряд в катушке). См. physics/losses.py.
+    r_eddy_coeff_ohm: float = 0.0
+
+
+def effective_resistance(inductance_model: CoilInductanceModel, params: StageCircuitParams, x) -> float:
+    """R контура = провод+ESR+ключ (r_total) плюс вихревые потери снаряда, когда он в катушке."""
+    return params.r_total_ohm + params.r_eddy_coeff_ohm * inductance_model.overlap_fraction(x)


 def derivatives(
@@ -37,10 +45,11 @@ def derivatives(
    q, current, x, v = state
    dlambda_di = inductance_model.dlambda_di(x, current)
    dlambda_dx = inductance_model.dlambda_dx(x, current)
+    r_eff = effective_resistance(inductance_model, params, x)

    v_c = q / params.capacitance_f
    d_q = -current
-    d_i = (v_c - current * params.r_total_ohm - dlambda_dx * v) / dlambda_di
+    d_i = (v_c - current * r_eff - 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]
--- a/src/gausse/physics/losses.py
+++ b/src/gausse/physics/losses.py
@@ -0,0 +1,55 @@
+"""Потери в железе снаряда: вихревые токи + гистерезис (первого порядка).
+
+Сплошной стальной снаряд в импульсном поле теряет энергию на вихревые токи
+(греется). Модель: снаряд — короткозамкнутый виток (вторичка трансформатора),
+резистивно-доминированный на частоте импульса. Отражённое в катушку
+сопротивление R_eddy = (ω·M)²/R_e добавляется в контур, ПОКА снаряд в
+катушке (× overlap(x)) — энергия честно уходит в нагрев снаряда, а не в
+кинетику. Это и есть тот канал потерь, из-за отсутствия которого КПД был
+нереально высоким.
+
+ВАЖНО (честно): это оценка первого порядка. Точные вихревые потери требуют
+решения уравнения диффузии поля в снаряде (скин-эффект, частотная
+зависимость). Здесь — сосредоточенный виток + характерная частота импульса
+ω=1/√(L·C). Геометрия вихревого контура (R_e) и коэффициент связи оценены
+приближённо; величина может отличаться в разы. Коэффициенты вынесены, чтобы
+их можно было уточнить по реальным замерам.
+"""
+
+import math
+
+from gausse.physics.constants import MU_0
+
+# доля потока катушки, реально сцепленная с вихревым контуром снаряда (k<1).
+# Снаряд у́же катушки и захватывает не весь поток — грубая оценка.
+EDDY_COUPLING_FACTOR = 0.5
+# гистерезис: энергия на площадь петли B-H за перемагничивание, Дж/м³/цикл.
+# Для конструкционной стали ~сотни Дж/м³; берём ~300 (справочная оценка).
+HYSTERESIS_LOSS_J_PER_M3 = 300.0
+
+
+def eddy_reflected_resistance_ohm(
+    slug_radius_m: float,
+    slug_length_m: float,
+    slug_resistivity_ohm_m: float,
+    mu_eff: float,
+    total_turns: int,
+    coil_length_m: float,
+    char_omega_rad_s: float,
+) -> float:
+    """Отражённое в катушку сопротивление вихревого контура снаряда (Ом)."""
+    a_slug = math.pi * slug_radius_m**2
+    turns_per_m = total_turns / coil_length_m
+    # взаимная индуктивность катушка↔снаряд (снаряд как 1 виток)
+    mutual = EDDY_COUPLING_FACTOR * MU_0 * mu_eff * turns_per_m * a_slug
+    # сопротивление сосредоточенного вихревого контура (кольцевой путь в снаряде)
+    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
+
+
+def hysteresis_energy_j(slug_radius_m: float, slug_length_m: float, n_cycles: float = 1.0) -> float:
+    """Энергия гистерезиса за выстрел (Дж) — вторичный канал, обычно << вихревых."""
+    volume = math.pi * slug_radius_m**2 * slug_length_m
+    return HYSTERESIS_LOSS_J_PER_M3 * volume * n_cycles
--- a/src/gausse/sim/stage.py
+++ b/src/gausse/sim/stage.py
@@ -30,6 +30,7 @@ from gausse.physics.inductance import (
 )
 from gausse.physics.constants import SWITCH_SURGE_FACTOR
 from gausse.physics.force import saturation_scale, solenoid_field_estimate_tesla
+from gausse.physics.losses import eddy_reflected_resistance_ohm
 from gausse.physics.sensors import (
    inductive_signal_peak_v,
    make_inductive_sensor_event,
@@ -139,6 +140,20 @@ def run_stage(

    demag = demagnetizing_factor_prolate(projectile.aspect_ratio)
    mu_eff = effective_permeability(projectile.material.mu_r, demag)
+
+    # вихревые потери в снаряде: отражённое сопротивление на характерной
+    # частоте импульса ω=1/√(L·C), действует пока снаряд в катушке (см. losses.py)
+    char_omega = 1.0 / math.sqrt(l_air_h * capacitance_f)
+    r_eddy_coeff = eddy_reflected_resistance_ohm(
+        slug_radius_m=projectile.diameter_m / 2,
+        slug_length_m=projectile.length_m,
+        slug_resistivity_ohm_m=projectile.material.resistivity_ohm_m,
+        mu_eff=mu_eff,
+        total_turns=geometry.total_turns,
+        coil_length_m=geometry.coil_length_m,
+        char_omega_rad_s=char_omega,
+    )
+
    inductance_model = CoilInductanceModel(
        l_air_h=l_air_h,
        coil_length_m=geometry.coil_length_m,
@@ -196,6 +211,7 @@ def run_stage(
        capacitance_f=capacitance_f,
        r_total_ohm=r_total_ohm,
        mass_kg=mass_kg,
+        r_eddy_coeff_ohm=r_eddy_coeff,
    )
    q0 = capacitance_f * stage.charge_voltage_v
    discharge_sol = solve_ivp(
@@ -219,8 +235,11 @@ def run_stage(

    energy_in_j = 0.5 * capacitance_f * stage.charge_voltage_v**2
    energy_remaining_cap_j = q_final**2 / (2 * capacitance_f)
+    # потери = I²·R_eff(x), где R_eff включает вихревые потери снаряда (× overlap)
+    overlap_during = inductance_model.overlap_fraction(discharge_sol.y[2])
+    r_eff_during = r_total_ohm + r_eddy_coeff * overlap_during
    energy_dissipated_j = float(
-        np.trapezoid(discharge_sol.y[1] ** 2 * r_total_ohm, discharge_sol.t)
+        np.trapezoid(discharge_sol.y[1] ** 2 * r_eff_during, discharge_sol.t)
    )
    kinetic_before_j = 0.5 * mass_kg * v_fire_mps**2
    kinetic_after_j = 0.5 * mass_kg * v_final**2
--- a/tests/test_stage_energy_conservation.py
+++ b/tests/test_stage_energy_conservation.py
@@ -38,7 +38,7 @@ def _switch(on_resistance_ohm: float = 0.02) -> SwitchSpec:
    )


-def _run(wire, switch, capacitor=CAPACITOR):
+def _run(wire, switch, capacitor=CAPACITOR, projectile=PROJECTILE):
    # turns_per_layer выбран так, чтобы длина катушки (~1.7см) была сравнима
    # со снарядом (2см) — иначе снаряд, войдя глубоко внутрь длинной катушки,
    # оказывается в плоской зоне перекрытия (dL/dx=0) и сила не действует.
@@ -53,7 +53,7 @@ def _run(wire, switch, capacitor=CAPACITOR):
        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)
+    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():
@@ -66,8 +66,16 @@ def test_energy_conserved_exactly_when_lossless():
        part_number="c-lossless", capacitance_uf=100.0, voltage_v=400.0, esr_ohm=0.0,
        max_current_a=300.0, price=300.0, source="test",
    )
+    # снаряд без вихревых потерь: огромное удельное сопротивление -> R_eddy≈0,
+    # чтобы «без потерь» действительно означало отсутствие всех каналов диссипации
+    no_eddy_steel = ProjectileMaterialSpec(
+        name="steel-no-eddy", density_kg_m3=7850.0, mu_r=200.0, b_sat_tesla=1.8,
+        price_per_kg=100.0, source="test", resistivity_ohm_m=1e12,
+    )
+    no_eddy_projectile = ProjectileConfig(material=no_eddy_steel, diameter_m=0.008, length_m=0.02)
    result = _run(
-        _wire(resistivity_ohm_m=0.0), _switch(on_resistance_ohm=0.0), capacitor=lossless_capacitor
+        _wire(resistivity_ohm_m=0.0), _switch(on_resistance_ohm=0.0),
+        capacitor=lossless_capacitor, projectile=no_eddy_projectile,
    )
    assert result.feasible, result.reason
    assert result.energy_dissipated_j == pytest.approx(0.0, abs=1e-9)