diff --git a/src/gausse/components/data/projectile_materials.json b/src/gausse/components/data/projectile_materials.json index b34542d..5d06c1b 100644 --- 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 Ом·м (чище -> ниже)" } ] diff --git a/src/gausse/components/schema.py b/src/gausse/components/schema.py index 615eae8..59746cc 100644 --- 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 diff --git a/src/gausse/gpu/batch_integrator.py b/src/gausse/gpu/batch_integrator.py index 1f54418..bbf5b60 100644 --- 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 diff --git a/src/gausse/physics/circuit.py b/src/gausse/physics/circuit.py index 61701df..83eb76c 100644 --- 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] diff --git a/src/gausse/physics/losses.py b/src/gausse/physics/losses.py new file mode 100644 index 0000000..335115f --- /dev/null +++ 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 diff --git a/src/gausse/sim/stage.py b/src/gausse/sim/stage.py index fbdccba..4e14934 100644 --- 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 diff --git a/tests/test_stage_energy_conservation.py b/tests/test_stage_energy_conservation.py index 0c83741..4c1b4c8 100644 --- 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)